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ABSTBACT 


This work describes the development of a computer based 
model that would allow for determination of the transient 
response characteristics of jet vanes of any size. The model 
used a thermal lump approach method, considering only 
conductive and convective heat transfer properties. A 
constrained optimizer was used to adjust the unknown variables 
until an adequate match was achieved between the calculated 
values of the energy balauice equations and the experimental 
data obtained from test firings of a rocket motor. The full 
scale modeling results were compared to previous quarter scale 
models in an atten^t to determine the applicability of the 
quarter scale results to full scale vanes. It was determined 
that the quarter scale model did not provide an accurate 
representation of the heat transfer process in larger scale 
vanes, although the full scale model provided a sufficient 
representation of the thermal response of the jet vane. 
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I. INTRODUCTION 


The work reported in this document is further research in 
the study of Thrust Vector Control (TVC) system used in 
missile trajectory guidance during initial phases of a 
missiles flight [Ref. 1]. TVC utilizes a system which 
deflects exhaust thrust to control the missile during launch 
and initial flight stages. For this phase of flight the vane 
is actuated from a stowed position and inserted into the 
exhaust of the rocket nozzle to provide the required control. 
Figure 1.1 illustrates the stowed and active positions of the 
TVC vane. 



Figure 1.1 STARS TVC Sto«#age Concept 
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The vane operates in an extremely harsh thermal 
environment. In order to use adequate materials for vane 
design the heat transfer characteristics of a selected 
material operating in such an environment needs to be 
resolved. In order to determine the best performance 
characteristics of different materials a mathematical model 
needs to be constructed. The model needs to contain both the 
known physical properties of the material and the unknown 
properties to be determined. It will utilize the experimental 
temperature time histories from rocket firings and adjust the 
unkno%ms physical properties until the predicted temperature 
time history and the experimental data are close in a least 
squares sense. 

Past work in the system identification process has been 
conducted using a quarter scale model to determine the 
convective heat transfer coefficients and to verify the 
applicability of the quarter scale results in predicting the 
heat transfer characteristics of a full scale protoype. The 
idea behind this report is to enlarge the model to full scale 
in order to obtain the heat transfer properties and to compare 
these results with those obtained from the quarter scale 
modeling process. 
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II. BACKGROUND 


A. LUMPED PARAMETER MODELING 
1. Vane Mbdeling Concept 

The jet vane that is part of the TVC system concept 
will be used as a deflection device to direct the thrust of 
the rocket exhaust. The vane will be encountering effects 
from such things as con^lex thermal changes, supersonic 
particle laden stream flow, shock waves impinging on the vane 
and phase variations resulting from vane ablation [Ref. 2:p. 
5]. The idea of vane modeling is a process by which some or 
all of these processes may be modeled in order to find the 
material that provides the durability and the necessary 
performance characteristics. Due to the mathematical 
complexity and the wide number of variables that may effect 
vane performance it is necessary to limit the model to a 
smaller scale. The main emphasis of this report is to limit 
the model to determining the thermal convective 
characteristics. This modeling concept involves development 
of a thermal model taking into account the conductive and 
convective effects of the material selected. Radiation 
effects will not be considered in the development of this 
model. The model uses the concept that the transient thermal 
response of the vane can be developed into a set of equations 
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that describe the heat transfer process occurring over a 
distinct period of time. 

2. Lumped Parameter Approach 

The thermal modeling approach taken will be to break 
the vane down into distinct thermal lunps (Ref. 2;p. 101. 
This approach lumps the heat transfer parameters into distinct 
nodes. The geometric data and physical properties of the vane 
material at are then used to determine the conductive 
resistance and capacitances at each node. The convective 
resistances between the exhaust plume and the material are 
essentially unknown and will be determined using the modeling 
process. The lumped parameter approach is essentially a 
division of the vane into critical locations< usually chosen 
to correspond to a thermocouple attachment point on the 
experimental model< in order to simply the modeling process 
(Ref. 2:p. 12]. This simplification enables the formulation 
of energy balance equations* mentioned previously* that 
describe the transient thermal response of the individual 
nodes in the vane (Ref. l:p. 4]. The energy balance equations 
take into account thermal energy .ito and out of the vane 
nodes and the temperature of the vane accounts for the energy 
stored at the nodes. 

3. Parametric System Identification 

The thermal model will be computer-based FORTRAN code 
using Parametric System Identification (PSD (Ref. 3:p. 2]. 
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PSI is a procedure in which the model parameters generate a 
temperature time history that predicts the experimental 
tenperature time history [Ref. 2:p. 6] . The process involves 
adjustment of the unknown parameters until a desired agreement 
is reached between experimental and model temperature 
profiles. The PSI process involves adjustment of the unknown 
values for the convective coefficients until the model 
tenperature time histories of selected nodes are matched with 
the experimental temperature time histories of the 
corresponding nodes on the actual vane. The PSI process 
brings together the mathematical model and the experimental 
data. The accuracy of the match between model and 
experimental data is limited to the complexity of the energy 
equations used to describe the model itself. The process can 
be a self adjusting one in that the model can be used to 
design the vane for experimental testing and then the 
resulting experimental data can be used to update the original 
model (Ref. 2:p. 6]. 

B. B^slc mosLim processes 

1. Ecruatioo Formulation 

The modeling of the TVC vane involves generating the 
energy balance equations for the lumped parameter model. The 
energy equations are formulated by considering the flow of 
energy into and out of a node. Figure 2.1 illustrates the 
energy flow concept of a vane node. 
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Figure 2.1 I3iezinal Node Energy Flow 
The energy balance formulation for Figure 2.1 can be 
constructed by considering the flow of energy into and out of 
the node. This energy balance is represented by, 


' CxRx C^Ro 


( 1 ) 


where the value of Tj>Ti>To. Starting with the determination 
of the resistances and capacitances at each individual node 
from the geometric data and the physical properties of the 
material is required. Prior to understanding the effects of 
sub-scale model the conductive and convective parameters first 
need to be addressed. The conductive resistance is found by, 

where L is the conductive path length, \ is the cross 
sectional area of the conductive path and K is the thermal 


6 






conductivity [Ref. 3:p. 7]. Likewise, the convective 

resistance can be found ly. 






{3) 

where A, 

is the convective surface 

area 

and h is 

the 

convective 

heat transfer coefficient 

[Ref. 

3:p. 7]. 

The 


thermal capacitance is found by, 

(4) 

where q is the material density, Cp is the specific heat and 
V is the volume [Ref. 3:p. 7]. It can be seen that these 
values are dependent on not only the material in which the 
vane is constructed but also the actual size of the vane being 
modelled. All geometric data is given in full scale and the 
physical properties are assumed constant no matter what size 
model is being considered. 

The vane model is generated by formulating the energy 
balance equations for each individual node of the vane. The 
equations will be dependent on the number of nodes used to 
model the heat transfer process occurring in the vane. To 
develop such a set of equations the nodes locations need to be 
established on the vane. The nodes chosen should correspond 
to a center of mass or the major sections of the object being 
modeled, in this modeling process the vane is broken down 
into four distinct sections. Figure 2.2 illustrates a vane 
with a four node configuration (Ref. l;p. 6). 
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Figure 2.2 Four Node Vane Model 
This four node select:ion is by no means a standard which must 
be adhered to. Multiple nodes may be selected depending on 
the complexity that is desired. 

Xn Figure 2.2 the vane is broken down into four distinct 
nodes. Node one is for modeling the tip, node two for the 
fin, node three for the shaft and node four for modeling the 
mount (Ref. l:p. 5]. The energy balance equations are 
generated by considering the heat transfer into and out of 
each individual node. Nodes one and two, two and three, three 
and four, are connected by conductive paths with a conductive 
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resistance linking the nodes [Ref. 3:p. 4]. The designation 
Rij is used to denote the conductive resistance between nodes 
i and j. Nodes one and two are linked convectively to two 
driving temperatures TRl and TR2. These are known as the 
recovury temperatures and are the source of the heating 
process [Ref. 3:p. 4]. TRl corresponds to the stagnation flow 
while 'rR2 represents the tenperature of the free stream flow. 
Tue convective resistances are represent by the designation 
Rpi, where the i represents the node that the convective 
temperature is driving. Each node contains a thermal 
capacitance represented by Ci where i is the node 
representation. The generation of the energy balance 
equations involves looking at each individual node and 
applying the law of conservation of energy [Ref. 3:p. 4). 
The conservation of energy for each node involves energy 
transfer in and out of the node. That is the rate of change 
of energy at a node is equal to the er.arny into the node minus 
the energy out of the node. The remaining energy at the node 
is what causes the temperature change at that particular 
location. The node is a point loca^-ion and the transient 
thermal heating is being considered. 

/^splication of this idea leads to the following equations 
[Ref. 3:p. 4); 







(5) 
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(6) 




In order to make conputer coding easier a sin^lification of 
the above equations is desired. The following characteristic 
rate coefficients are defined (Ref. 3:p. 9]: 



1 



(9) 

«3a“ 

1 



(10) 

^40“ 

1 

Jx. . 

CxRrt 

b,,- ^ 

(11) 


A further simplification can be made by combining coefficients 
for the same temperature at a particular node. This 
simplification allows each node temperature to have a 
single coefficient making it possible to put the energy 
balance equations into matrix form. This is done as follows 
(Ref. 3;p. 5]; The energy balance equations now become (Ref. 
3 s p • 5 ] 
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^33 “^32 ■*■^34 ( 13 ) 

^2®^21^l“*22^2‘'‘®23^3''’-^2^Jl2 

(16) 

^ 4=^43 ^ 3 - 344^4 ( 17 ) 


These equations may be put in matrix form to give a clearer 
understanding of there use in the computer coding 
[Ref. 3:p, 5], 

—Qn <*12 0 0 ' " Ti ' ' frii 0 0 0"’ 7 /ji 

ti JJ, <*ia ”<*22 <*23 0 ^2 4 . 0 622 0 0 Tr 2 

T3 0 03J —033 ^34 Ta 0 0 0 0 0 

ti, L 0 0 043 -<*44 J I T4 . . 0 0 0 0 . . 0 . 

This is a compact notation commonly known as a state space 
equation and may be written as, 

jic°AJC'^Bu ( 19 ) 

Equations of this form are common and a wide variety of 
package are available for solving linear differential 
equations. 


11 





2. Previous Slodels 


Work on the jet vane thermal model was started by Nunn 
and Kelleher in 1986 at the Naval Postgraduate School [Ref. 
5] . Further research into the model was conducted by Nunn 
[Ref. 2], Hatzenbuehler [Ref. 4] , Reno [Ref. 1] and Driels 
[Ref. 3]. Prior work has concentrated on development of a 
quarter scale model and its applicability to full scale vanes. 
Reno used both 3 and 4 node approaches in modeling the heat 
transfer process of the vane. The experimental data used for 
the comparison process, on the quarter scale model, was from 
the shaft and mount locations [Ref. l:p. 7]. The points were 
chosen were due to the limited size of the experimental veuie 
for thermocouple attachment. Those points being the back of 
the shaft and mount. The 4 node vane model in figure 2.1 
illustrates the set-up used by Reno along with the energy 
balance equations. The unknovm variables that needed to be 
determined were Rp 2 # ^ 34 * K 4 Q and C 4 [Ref. l:p. 10] . In the 
simplified energy balance equations these coincide to a 34 , a 43 , 
a 4 Q and b 22 < The geometric data used for the full scale vane 
is given in Table 1 [Ref. 3:p. 8 ]. 

TABLE 1 



tip to vane to shaft 

V, cm^ 

A^, cm^ 

L, cm 

2.6 52 23 

5.9 5.2 

5,0 6.0 
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The physical data for the vane, which was a tungsten matrix 
impregnated with copper, is p = 18310 kg/m^ K = 173 W/m-K, 
and Cp = 146 J/kg*K [Ref. 3:p. 7]. The recovery tenperatures 
used were TRl = 2360 K and TR2 = 2260 K. The input thrust 
profile used was simulated as a ramp function. The recovery 
temperatures, TRl and TR2, contained in the input vector u, of 
the matrix form, drive the system, and are the product of the 
steady state values and the thrust function. Figure 2.3 
illustrates this functional input [Ref. 3]. 



Figure 2.3 Rasp Function for Recovery Tenperatures 
Reno utilized a personal computer with a program System Build 
and System Identification developed by Integrated Systems, 
Inc. [Ref. l:p. 10] . The identification process used MAXLIKE 
which is a function of MATRIX X [Ref. l:p. 10]. This program 
allowed the four unknown parameters to vary until the 
temperature profiles generated matched those obtained from the 
experimental firings. The following results were obtained: 
bjj o 0.2718 s’‘, aj 4 = 0.299 s'‘, a*, = 0.1322 s‘‘ and 
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340 = 0.1018 s'^ [Ref. I:p.l0]. These parameters correspond to 
the values of [Ref. l;p. 11]: R„ = 1.67 K/W, R 34 = 3.33 K/w, 
R 40 = 4.33 K/W, C 4 = 2.27 J/K. There was a good match between 
the generated and experimental temperature time histories of 
the mount and the shaft. Figure 2.4 demonstrates the close 
proximity between the shaft and mount profiles (Ref. l:p. 12] . 



TtNF. (K) 

Figure 2.4 Quarter Scale Tenperature Profiles 
Determination of the unknown paramet< s is based upon a 
temperature profile match between the mount and shaft. The 
temperature profiles ten^eraturcs and Tj are the 

mathematical equation profiles and have no experimental data 
to be matched to in the quarter scale modeling process. 
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Work done by both Reno (Ref. 1] and Driels [Ref. 3] 
indicated that although the convective heat transfer 
coefficient at the tip node was an unknown variable, in the 
quarter scale modeling process, it had minimal effect on the 
outcome of the previously identified unknown values and thus 
bii was held constant. Table 2 illustrates the variation in 
the four unknowns for various values of bn [Ref. 3:p. 12]. 

TABLE 2 


:l9i 

a34 

a«3 

^40 

^22 1 

5.35 

0.483 

0.152 

-0.167 

0.143 1 

10.0 

0.447 

0.152 

-0.168 

0.162 

20.0 

0.44 

0.152 

-0.169 

0.174 

1.0 

0.51 

0.153 

-0.164 

0.056 

0.5 

0.47 

0.153 

-0.163 

0.019 


The results reported in Table 2 are those obtained by Driels 
[Ref. 3] and do not exactly match the results obtained by Reno 
[Ref. 1]. The tert®)erature profiles obtained by Driels 
coincide with those displayed in Figure 2.4. 


C. MODEL SCALING 

1. Purpose of Quarter Scale Model 

The quarter scale model was developed to predict the 
behavior of the heat transfer pareimeters of a full scale 
prototype. The quarter scale model quantitatively provides 
for testing on a smaller scale which in the long term is more 
cost effective than constructing one of full scale. This 
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savings come from the reduced size of the rocket motor, the 
reduced use of data acquisition equipment, environmental 
requirements that have to be met when firing full scale that 
do not apply to the quarter scale, manufacturing of quarter 
scale vanes, material requirements, etc. The quarter scale 
modeling process, if adequate, could provide a quick, simple 
form of modeling for a larger vane without the added 
difficulties. 

2. Quarter to Full Scale Model Predictions 

It has been demonstrated that a sub-scale model can be 
generated and subsequent experimental data obtained in order 
to predict the heat transfer parameters of a given material. 
A set of convective heat transfer coefficients can be 
determined from the results of such models. The accuracy of 
these results and there applicability to a full scale vane has 
yet to be determined. Direct scaling of the geometric 
properties is achieved and the physical properties such as 
density, thermal conductivity and specific are assumed 
constant for both sub-scale and full scale vanes. The 
difficulty arises in the effects of the jet plume on the 
different sizes of vane. As mentioned earlier the conductive 
properties rely on both physical and geometric data. 
Variation in the geometric size would have a profound impact 
on the transient heat transfer process occurring not only 
inside the vane but due to the larger surface area the effects 
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of the convective heat transfer process. A direct scaling up 
of the quarter scale model to full scale may be inadequate to 
properly model the effects of a large vane. Increases in the 
rocket chamber temperature thus an increased recovery 
tenperature may also influence the energy transfer process. 
The purpose of this report is to develop a full scale model, 
compare the results to those obtained by the quarter scale and 
determine if the results obtained from the sub-scale tests 
reflect the actual heat transfer process of a full scale vane. 

3. Scaling Analysis 

To scale the resistances and capacitances a reduction 
in the linear dimensions of the vane must be calculated in 
order that the model reflects the actual size of the 
manufactured vane in the experimental phase [Ref. 3:p. 8], A 
scale factor SF, where SF<1, is introduced to reduce the 
geometric dimensions of the vane and subsequently altering the 
resistances and capacitances [Ref. 3:p. 8]. The scale values 
can be determined by manipulating equations (2), (3) and (4). 
When manipulating these equations the physical properties such 
as the thermal conductance k, the specific heat Cp, and the 
density q remain constant and do not change. The equation 
manipulation results in, 


' KA, 




L-SF ^ 


KA/SF^ 




_L _ 3 ^ 

KA,SF 


( 20 ) 
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( 21 ) 



where Rf is the full scale resistance and R, is the sub-scale 
resistance. An important point to note is that the scale 
factor SF being less than one causes the resistance to 
increase during the scaling process. An exanple of this is 
the quarter scale models previously constructed where the 
scaled resistance R, is equal to the full scale divided by 
0.25 (1/4 scale modeling). This results in the scaled 
resistance being four times the full scale value. The thermal 
capacitance has the same linear dimension scaling as the 
resistance did, which is given by, 

( 22 ) 

C^*>Cf’SF^ ( 23 ) 

where C( is the full scale capacitance and as indicated is 
the sub-scale capacitance (Ref. 3:p. 8). There is no direct 
scaling of the convective resistance in the model due to the 
fact that the heat transfer coefficient h is unknown. The 
only scaling used on the convective values is if convergence 
has been achieved and the value of the convective resistance 
is needed. The scaling process for the convective resistance 
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is similar to those previously stated and is given by [Ref. 
3:p. 81, 



Although there is an increase in the resistance due to scaling 
and a decrease in the capacitance, the CCTtbined effect of 
these two values must be emphasized. The product of RC is a 
time constant for a particular driving temperature [Ref. l:p. 
7]. Scaling of the values of R and C results in a change in 
the time constants for the energy balance equations from the 
full scale to the sub-scale. The resultant change is 
demonstrated in the following formula; 

- R,C,°RtCf'SF^ (25) 
SF 

The time constant for the sub-scale values decreases as the 
product of the full scale time constant times the square of 
the scale factor. It will be seen that in formulating the 
energy balance equations for the mathematical model that the 
scaling of the time constants will have a profound effect on 
the conductive and convective values. 
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III. FULL SCALE M}DELING 


A. SIX NODS MODEL APPROACH 

1. Concept 

As discussed previously, prior quarter scale modeling 
approaches involved using three and four node models. 
Placement of the node points in the model were dependent on 
the center of mass at determined critical locations and the 
availability of experimental data from a quarter scale test 
firing. The experimental data coming from the shaft and mount 
locations. These temperature profiles were matched by the 
matheinatical model and the results used in determining the 
convective heat transfer coefficient for the fin section, the 
unknown resistance between the mount and shaft, and the 
unknown capacitance of the mount. A single value of 
convective heat transfer coefficient was obtained from the 
quarter scale tests but due to the lack of experimental 
temperature profiles from that particular location the 
applicability of this value to a larger vane was in question. 

The basis for increasing the nuitdaer of nodes when 
attenpting to model the full scale vane was two fold; 1) more 
node locations would give a more accurate representation of 
the actual heat transfer process going on in the vane and 2) 
these added node locations had experimental data available 
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which allowed for a more diverse set of temperature time 
profiles that could be used in constructing the mathematical 
model. The 21 locations that experimental data was taken from 
allows a greater flexibility in the selection of temperature 
profiles to be matched. Although Nunn [Ref. 2:p 12] 
recommended that for simplicity and efficiency the minimum 
number of nodes should be used, to accurately model the heat 
transfer process of a full scale vane a larger number of nodes 
was needed. The modification was made by adding two more node 
locations to the four node model in the fin area. Figure 3.1 
shows the six node model design. 
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Figure 3.1 Six Kode Vane Model 
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2. Geometric Modification 

The addition of two nodes to the fin in the six node 
model require a modification of the geometric properties in 
order that the proper volumes, areas and lengths are used in 
calculating the resistances and capacitances. These new 
values are given in Table 3, 

TABLE 3 



tip to node2 to node3 to node4 to shaft 

V, cm^ 

A*, cm^ 

L, cm 

2.6 12.86 17.15 21.44 23.0 

5.6 5.9 6.2 5.2 

2.774 2.774 2.774 5.0 


The values are based on those obtained from Table 1 [Ref. 3:p. 
8] and the four node model design. The parameters for the 
node volumes were obtained by taking the vane volume and 
estimating a 20% change in the volume from the rear of the 
vane, node four section, to the forward part of the vane, node 
two section, with the rear being the largest and the basis for 
the calculations. The cross sectional areas were determined 
in the same manner with the area at node two being considered 
to be 20% less than node three and node four cross sectional 
area to be 20% greater than that of node three. The cross 
sectional area at node three was known because it corresponds 
to node two of the quarter scale four node model. The quarter 
scale model parameters were given in full scale terms and 
scaled down in the program formulation. 
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3. Energy Balance Equation Foinnulation 

Formulation of the energy balance equations is done by 
applying the law of conservation of energy to the six nodes of 
the vane. Application of this law leads to the following, 
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As discussed in Chapter II this series of equations can be 
simplified by defining a set of characteristic rate 
coefficients. These coefficients relate the properties of the 
vane to the time constants. The simplification involves both 
a conductive coefficient labeled a^^j and a convective 
coefficient labeled b^j. This simplification leads to, 
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These simplifications now allow the energy balance equations 
to be put into matrix form. The modification provides an 
easier approach to the program formulation. The matrix 
equation description was described in Chapter II and with the 
application of the characteristic rate coefficients the 
lumping of the common ten^erature coefficients in each energy 
balance formulation the following equations can be used to 
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generate the FORTRAN code used to solve these differential 
equations. The equations become. 
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4. Program Formulation 

The energy balance equations are a set of liner 
ordinary differential equations with some of the rate 
coefficients as unknowns. The process to evaluate these 
unknowns is to start with a set of initial conditions for 
these unknowns, solve the differential equations over a period 
of time to obtain a tenperature time history for each node. 
We then con^are the experimental data from a set of 
corresponding nodes and the results obtain from the model 
progreun and calculate the difference of the sum of the squares 
between the two temperature profiles. The programs optimizer 
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adjusts the unknovm variables and solves the differential 
equations again to obtain a new temperature time history. 
This process is illustrated in Figure 3.2 



Figure 3.2 Block Diagram of Modeling Program 
The blocks labeled Pre.dat, Trans.for, TV.mat,TM.mat and 
Pro.m will be discussed in the following section. The program 
PIO reads the experimental data from the Datam.dat file. It 
then reads in the full scale geometric properties of the vane 
and the physical properties of the material the vane is 
constructed of from Xnput.dat. It calculates the resistance 
and capacitances for the nodes and the rate coefficients a^j. 
The initial conditions are set for the unknown parameters and 
PID then calls the optimizer subroutine. The optimizer 
subroutine intern calls another subroutine called TEMP which 
calls the needed equation solver and function program that 
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provides the differential equations. DIVPRK (Double-precision 
Initial Value Problem Runge-Kutta) is an IMSL routine used to 
solve linear ordinary differential equations with the given 
set of initial conditions. The subroutine that supplies the 
differential equations is called FCN. The values of the 
driving temperatures are supplied by the main program PID and 
are common to all the subroutines. The results of DIVPRK are 
returned to the subroutine TEMP as a temperature time history. 
TEMP then calculates a sum of the squares error between the 
experimental data read into the program PID and the results of 
the equation solver. This error is returned to the optimizer 
program where the optimizer adjusts the unknown parameters and 
repeats the process until a set of convergence criteria is 
met. 

Once the program converges the results are vfritten to two 
files Result,dat and Temp.mat. Results.dat contains the 
calculated values of the convective resistance and the 
convective heat transfer coefficient. The temp.mat file 
contains the experimental and model temperature time 
histories. This file can be read into MATLAB to be graphed, 
which gives a visual representation of how close the match is 
between the two profiles. 

5. Experimental Data Processing 

The experimental data used in the identification 
process is obtained from the Naval Weapons Center (NWC), China 
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Lake, California. data is received on a o 'ter readout 
sheet for the different thermocouple locations on the 
e*<perimental vane. The temperature is measured in degrees 
Fahrenheit. The computer program needs the temperatures in 
degrees kelvin and referenced to the ambient or reference 
temperature. This referencing of the temperature profiles 
enables the initial conditions of the equation solver to be 
set to zero. In order to use the experimental data it must 
first be preprocessed. The preprocessing of the data is 
illustrated in the block diagram in Figure 3.3. 



Figure 3.3 Experimental Data Preprocessing Diagram 
The Pre.dat file is first produced by manually selecting the 
series of tenperatures of a specified time increment. The 
separation of the temperatures is up to the individual. The 
time interval selected when recording the experimental 
temperatures must also be placed into the TEMP subprogram for 
the differential Equation solver DIVPRK. This time interval 
allows the eciuation solver to calculate the temperature time 
histories over the S6une time increments as the experimental 
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data being used. Failure to match these time increments may 
cause erroneous values when trying to match experimental and 
calculated values. The values are entered into the file 
Pre.dat manually as a nxl matrix. The raw data from the 
Pre.dat file is read into a preprocessing program called 
Trans.for. This smoothing program involves adding the two 
ten^eratures prior to the one being averaged, the temperature 
being averaged and the two temperatures subsequent to that 
particular tenperature. Since there are no temperature values 
prior to the first two and subsequent to the last two, these 
values are taken as is. For 61 temperature values over a 
given period of time temperatures one and two, and 60 and 61 
are taken as is. The temperatures 3 through 59 are averaged 
values once the Trans.for program has been run. 

The number of files produced from Trans.for may vary 
depending on the number of experimental temperature profiles 
used. The basic files produced are Datam.dat and TV.mat. The 
IM.mat file is only used if experimental data is being 
utilized from the mount. The files TV.Mat and IM.mat are read 
into a MATLAB program called Pro.m. This program normalizes 
the tenqperatures to the ambient and converts the temperature 
from degrees fahrenheit to degress Kelvin. The output of this 
program is a file called Datam.dat which is a nxl matrix. 
Datam.dat is the smoothed file containing the time* temperature 
profiles of the experimental data. At this point the 
Datam.dat file is ready to be used by the main program. 
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6. Forward Model Matching 

a. Coapsir±aon of Steady State T&igperatiires 

A forward model is the reverse of the vane 
optimization model. It takes a set of values and generates a 
set of time-temperature profiles for a given set of data. 
There is no optimizing of the variables. All the rate 
coefficients are constants and the profiles generated are 
based on the solution to the differential equations using 
these constants over a given period of time. A forward model 
was constructed from the program PID to check the types of 
time-tenperature profiles that would be generated by these 
equations given a set of values for the rate coefficients. To 
construct the forward model from the main program PID the 
optimizer subroutine was eliminated and the PID program called 
the TEMP subroutine directly. This forward model allows the 
user to solve a set of equations over a given period with a 
set of rate coefficients, and b^^, remaining constant. A 
set of time-tenperature profiles was constructed from the full 
scale energy balance equations using the final optimized 
values for the unknowns that Reno found from the quarter scale 
modeling process [Ref. l:p. 10]. The values were scaled up 
from quarter to full scale and inserted into the PID program 
for the unknown values at the corresponding locations. The 
values of b 22 and b 24 did not correspond to anything that was 
done by Reno (Ref. 1], so the values for the convective rate 
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coefficients were all set equal to each other. The forward 
model was run on the six node vane model for both quarter and 
full scale using the following values: Quarter Scale - 
a34=0.299, a43=0.1322, a4c=0.1018, bu=30.163, b22=0.271 Full 
Scale - a56=0.01869, a65=0.008262, aeG=0.006362, bii=1.885, b 22 = 
b 32 = b 42 =0.017. If the thrust profile is allowed to continue 
indefinitely the steady state values of tenperature should be 
the same for both the quarter and full scale values on the six 
node vane model. These values can be determined analytically 
by considering a voltage divider across the resistance path to 
ground. The convective temperature sources act as the supply 
voltages in such a case. The analytical calculations for both 
sets come up to the same steady state temperature values for 
both scaling sizes as shown in Figure 3.4. The main 
difference is that the quarter scale reached its steady state 
ten^erature at a much quicker rate than the full scale. The 
quarter scale reached steady state in approximately 20 seconds 
where as the full scale had not reached its steady state 
ten:«)eratures even after 160 seconds. Figure 3.4 shows the 
quarter scale response of the node temperatures over a 160 
seconds firing period. 
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Figure 3.4 Quarter Scale Forward Model Tai^erature Response 
The response curves in Figure 3.4 are for a six node vane 
model with the nodes at the shaft and mount being neglected. 
The analytical values of the steady state temperatures at 
nodes one through four are; Tl=2340 K. T2a2210 K, T3al639 K 
and T4s926 K, obtained from the figure above. The values of 
the rate coefficients are set to those corresponding to the 
full scale vane and the forward model program is executed 
again. The time temperature profiles for the full scale are 
listed in Figure 3.5 over the same 160 second firing period. 
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Figure 3.S Full Scale Forward Model Tesoperature Response 
The analytical calculations yield the same steady state values 
for nodes one through four. The variation is in the amount of 
time needed to get to the final value. The full scale 
requires a much greater period of time to reach the steady 
state than the quarter scale does, although the full scale 
will reach the same steady state temperature. This goes back 
to the time constants of the two different scale models. The 
rate characteristics decrease for the full scale model but 
these are the reciprocals of the time constant which actually 
increase for the full scale model. 
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b. Rasponse Time Characteristics 

A check needs to be made to see if the response 
characteristics of the energy balance equations for the six 
node model are similar to those of the four node model 
obtained from Reno's work [Ref. 1]. The forward model will be 
utilized with the values of the rate coefficients set equal to 
those in Reno's thesis [Ref. l:p. 10]. The thrust profile 
will be the same as that used by Reno and so will the recovery' 
temperatures, with Tri= 2360 and Tr 2 =2260. The time-temperature 
profiles obtained by Reno for a four node quarter scale model 
are shown in Figure 3.6. 
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Figure 3.6 Temperature Histories, Four Mode Model 
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The results in Figure 3.6 are based on an optimization program 
but if the energy balance equations are correct the final 
values used by Reno should be able to be inserted into the 
forward model generating a set of temperature profiles similar 
to those shown above. Figure 3.7 shows the results of this 
forward model using the six node energy equations with Reno's 
rate characteristic coefficients. 



TUIE (••€) 

Figttre 3.7 Six Node Quarter Scale Forward Model Response 
Comparison of the temperature profiles of Figures 3.6 and 3.7 
indicate they are very similar with the temperature profiles 
of T2 from the four node and T3 from the six node 
corresponding to each other. T3 from the four node 
corresponds to T5 from the six node and T4 from the four node 
corresponds to T6 of the six node model. 
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7. Model Simulation ZZSSQ Optimizer 

a. UnknocmB and ABBvapticma 

The first series of runs utilized an optimizer 
knovm as ZXSSQ which is an old IMSL routine that is stored in 
a separate library. It must be noted at this point that the 
thrust profile for the full scale modeling process is not a 
ramp function as those conducted in the quarter scale process 
but is being simulated as a step response. The thrust turns 
on at approximately time zero and shuts off approximately 6.7 
seconds later. Prior to starting the PSI process using the 
program formulation provided, the unknowns must be identified 
from the energy balance equations and a set of assumptions 
must be established. The unknown values that need to be 
determined are the convective resistances from the free stream 
to the vane, Rp 22 « Rp 23 ' % 34 , conductive resistance 
between nodes five and six, Rgg, the conductive resistance 
from node six to ground, R^q, and the capacitance of node six, 
Cg. Although the tip convective resistance and bj^j^ are 
unknown, the first series of the identification process 
assumed this value to be constant as it was in the quarter 
scale siimilation conducted by Driels [Ref. 3} and Reno (Ref. 
11 . 

Two simulations were conducted; the first with the 
assumption that the convective resistances to the free stream 
were dependent and equal, that is Rp 22 *% 23 "% 24 ' makes 


36 







b 22 =b 32 =b 42 * The second simulation made the assumption that the 
convective resistances were independent and not equal. 

Jb. Utilization of Experimantal Vane Tea^ Profiles 

Prior to model simulation the node locations for 
the use of experimental data needs to be selected. The 
quarter scale modeling done in previous work has utilized 
experimental ten^erature profiles from the shaft and mount 
locations. For the full scale vane esqjerimental data is 
available from the fin area of the vane. A set of simulations 
will be conducted using the same set of locations as that done 
by the quarter scale modeling process. Then the experimental 
data will be changed and the temperature profiles from nodes 
two, three and four will be utilized as the experimental time- 
ten^erature profiles to be matched. 

c. Results 

A series of simulations were conducted, with the 
initial conditions set at zero, varying the assumption of 
dependence and independence of the convective rate 
coefficients and varying the experimental data used in the 
identification process. All the simulations had the same 
outcome, there was no convergence. The values of as« and 
increase substantially in value causing a fatal error three in 
the differential equation solver giving the error message that 
the differential equation problems became too stiff. An 
attempt was made to alter the starting initial conditions in 
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the hope that the values could be placed closer to their final 
values. This was a manual attenpt to correct the sudden 
change in the values for the unknowns. Varying the initial 
conditions had no effect on the outcome. The problem still 
failed to converge. 

8. Model Simulation DBCLSF Optimizer 
a. UsB of Constrained Optimizer 

Since the values of some of the unknowns would jump 
substantially causing the program to fail, a constrained 
optimizer was used in an attempt to hold these values within 
a certain set of limits. The constrained optimizer chosen was 
DBCLSF which is a current IMSL subroutine. The reasoning 
behind this was that maybe the jump in the unknowns was caused 
by a local minima and when the optimizer, ZXSSQ, made a large 
adjustment it caused the problem to become to stiff. The 
constrained optimizer would keep the values of the unknowns 
low enough to prevent the differential equation solver from 
becoming to stiff and allow it to continue to solve the set of 
equations given. 

Jb. Boundary Limitations 

The constrained optimizer allows the user to select 
a variety of boundary conditions. This optimizer, DBCLSF, 
permits the selection of boundary limits that are all 
negative, all positive, or manually set either to each 
individual unknown or as a blanket setting for all the 
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unknowns. A simulation was conducted with the boundary limits 
set to all positive and then to a set with the lower boundary 
being set at 0.0001 and the upper boundary set at 100. 

c. Rasulta 

For the simulation, with the set of initial 
conditions at zero, there was no convergence with either the 
boundary limits set all positive or for the upper and lower 
boundary limits. Again the initial starting values of the 
unknowns were changed from zero to various values in an 
attempt to avoid any possible local minima. The results were 
the same as before the program failed to converge. The 
program had an arithmetic fault or just continued on 
indefinitely with the values of the unknowns reaching the 
boundary limits and remaining there. The program is trying to 
match a series of time-teirperature profiles from the fin area 
of the vane. The program does not converge due to the 
substantial jun«) in the unknowns at the mount location. From 
the experimental data the temperature profiles at the shaft 
and mount locations vary little from the ambient. Elimination 
of the mount node, node six, should have little effect on the 
predictions for nodes two, three and four. This coupled with 
the inconsistencies in use of mount configurations leads to 
the modification of the original six node model to a five node 
model. A model without a mount node and with the shaft node, 
node five, going directly to ground. 
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B. FIVE MODE MODEL APPROACH 

1. Concept Behind Mbdlfication From Six Node Model 

The computer simulations conducted on the six node 
model encountered the same problem in each simulation. The 
values of the unknowns ass and a^s, corresponding to the mount 
location, would juir^) a large amount in value causing the 
program differential equation solver to experience a fatal 
error, indicating that the equations were too stiff. Changing 
the relationship between the convective rate coefficients, bj^, 
had no effect on the resultant outcome and the jump in these 
unknowns prevented the convergence of the program. These 
unknowns contribute to the tenperature profiles primarily at 
the mount location. The convective heat transfer coefficients 
that were to be identified were between the free stream and 
nodes two, three and four of the vane and at this point we 
were not concerned with the mount. Examination of the 
experimental time tenperature profiles at the mount location, 
for the full scale vane, revealed that the temperature at 
these locations had little variation. The temperature varied 
less than 25°F over the entire >.7 second rocket firing period 
as conpared to the temperature variation at nodes two, three, 
and four which varied from 900 to 3000® F during the same 
rocket firing time. In contrast to the quarter scale modeling 
the temperature variation in the mount, node four, was 
significant. The size of the quarter scale vane limited the 
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locations for attachment of thermocouples to the shaft and the 
mount. This limited the experimental temperature data to only- 
two locations, the mount and shaft, making it necessary to use 
the mount temperature profiles in the matching process. The 
full scale vane is larger in size and provides greater space 
to place thermocouples and collect the experimental data. 

The model was constructed to determine the convective rate 
coefficients from the free stream to the vane. Due to the 
small variation of temperature in the mount and the enormous 
change of the unknowns using the mount node it was decided 
that a modification needed to be made to the six node model. 
This modification was the elimination of the mount node, which 
eliminated the two unknowns, age and a^s that were causing the 
program to fail. There was also another motivation for 
eliminating the mount node, which was the inconsistency in the 
use of the mount apparatus for all the experimental rocket 
firings (Ref. l:p. 13]. The modification to a five node model 
is illustrated in Figure 3.8. 
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Figure 3.8 Five Node Vane Model 
2. Energy Balance Equation Formulation 

The formulation of the energy balance equations for 
the five node model was done using the same process as that 
for the six node with only minor changes. The modification 
involves eliminating node six which in turn eliminates the 
last equation from the six node model energy balance. There 
was also only a minor change to the equation for node five of 
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the original six node equations. The following equations are 
applicable to the five node vane model; 


i*! = Ty + ai 2 ^2 Tgj 

(45) 

^2 ~^2l ~^22 ^2 ■‘■^23 ^3 ■*■'^2 

(46) 

^3 "^ 32^2 “^33 ^3 '*'^34 ^4 ■*‘■^2 

(47) 

^4 "^43 ^3 “^44 ^4 ■*‘345 ^5 ■'■^42 ^32 

(48) 

^ 5 “® 54 ^ 4 ~® 55^5 

(49) 

C^R^g ^ 55 “^ 54'‘‘350 

(50) 


Reduction of the model to five nodes eliminates the need to 
identify the unknown conductive resistances from the shaft to 
the mount, the mount to ground, Rjo, and the capacitance 
of the mount node, Cj. The experimental data for the five 
node full scale modeling process will be from the thermocouple 
attachment points that correspond to nodes two, three and four 
of the vane model. 

3. Program Modification 

The main program PID needs to be modified from the six 
node model to that of the five node. This is done primarily 
by making an alteration in the FCN subroutine which generates 
the differential equations to be used. The only other 
modification that needs to be made to the PID program is the 
set up of the size of the dimensions for the rate coefficients 
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and the setting of the initial conditions. Identification of 
the unknowns will be discussed in the following section. The 
FORTRAN coding used to generate the simulation program was 
constructed to be flexible in that it allows for changes in 
the vane nodes without major program modification. 

4. Model Simulation ZXSSQ Optimizer 
a. Uziknovns and Assvs^tions 

The modeling process for the five node model will 
be set up in the same way the six node model utilized the 
ZXSSQ optimizer. Prior to starting the simulations the 
unknowns for this model need to be identified. As before the 
unkiiowns are the convective rate coefficients, b22/ b32 and b42, 
along with the conductive rate coefficient, asQ. In the 
quarcer scale modeling it was demonstrated that the convective 
rate coefficient from the tip, b^, had little effect on the 
determination of the unknown parameters. Without knowing the 
effect on the full scale model the value of b^ will be 
considered a constant for the first series of simulations. 
The assumption of dependence between the convective rate 
coefficients was applied to the first series of simulations 
making b22=b}2=b42. A second set of simulations will be 
conducted making the values of the convective rate 
coefficients independent of each other and making the tip 
convective rate coefficient a variable. 
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b. Ea^riaental Data Used 

The experimental data used for the five node vane 
model will be from the thermocouple locations that correspond 
to nodes two, three, and four of the model. For the 

simulations involving the use of the ZXSSQ optimizer, an 
attempt will be made to match all three time-tenperature 
profiles. A point to be noted is that the values of the 
recovery temperatures that are to be utilized are different 
from those used in the quarter scale model due to the 
different rocket motor used in the test firings. The 
tenperatures used to drive the heat transfer process were 
Tri= 2970 and Tr 2=2870 Kelvin, respectively [Ref. 2:p. 44]. 
Referencing these tenperatures to the ambient, as was done 
with the experimental node tenperatures, would lead to 
recovery temperatures being normalized to values of Tri= 2670 
and Trj» 2570 Kelvin, 
o. Results 

As seen with the six node modelling process, using 
the ZXSSQ optimizer, the model failed to converge. Changing 
the assumption for the dependence of the convective rate 
coefficients, which made the values of the convective rate 
coefficients independent variables had no effect, and the 
program still failed to converge. An attenpt was made to 
match a single node time temperature profile but this also 
proved fruitless and again the program failed to converge. 
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From the foregoing results it is evident that the use of the 
ZXSSQ optimzer for the full scale modeling process is 
inadequate. There is no provision in the ZXSSQ optimizer to 
control the values or range of the unknown variables. 

5. Mbdel Simulation DBCLSF Optimizer 
a. Conatra.in3d Modal Simulation 

As with the six node model simulations, a 
constrained optimizer replaced the present optimizer, ZXSSQ. 
This was an attempt to hold the unknowns within certain limits 
permitting the optmizer a chance to avoid or recover from any 
possible local minimum and prevent the unknowns from becoming 
so large that the program fails to converge. The DBCLSF 
constrained optmizer utilized will be from the IMSL library. 
The first simulation assumed that the convective rate 
coefficients bjj, Djj and b„ are equal and that b^ was a 
constant and attempted to match the ten^erature profiles of 
nodes two, three, and four. This simulation also failed to 
converge. The assumption that bn was a constant was changed 
and bxi was allowed to vary. This simulation also failed to 
converge. 

b* Single Soda Temperature Profile Modaliag 

A different approach was taken to try and establish 
a five node working model of the vane. Instead of trying to 
match all three temperature profiles the model was altered to 
try and match only a single node time-teitq;>erature history then 
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work up to two nodes and then a full three node matching 
process. The assunptions, for the single node tenperature 
matching process, that any of the unknown variables were 
dependent were abandoned and each unknown variable was 
considered independent of the others. 

(1) Node Two Tenperature Modeling. The modeling 
simulations began by trying to match the time-tenperature 
history of only a single node, node two. This was attenpted 
in the hope that if only one profile could be matched that by 
examination of the profiles generated for the other nodes a 
slight adjustment could be made in order to obtain a match 
between all three node time tenperature profiles. The program 
was modified to match the profile of node two and the boundary 
limits were set to 0.0001 for the lower boundary and 100 for 
the upper boundary, with recovery temperatures of T,*j»2670 and 
Tm» 2570, and the initial conditions set to zero. The 
simulation converged with a sum of the squares error of 12.01, 
with values for the unknowns b„=3.0833, b,,*0.0255, b,a*0.0492, 
and b^,*0.2496, The value for a« reached the lower limit of 
0.0001 set in the constrained optimizer and did not vary from 
that point. Figure 3.9 illustrates the time-temperature 
profiles for this sisoilation. 
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Figure 3.9 Time-Tenverature Profiles Node Two Match 

(2) Node Three Temperature Modeling. The main 
program PID was modified to match experimental data for node 
three only. The program simulation converged with a sum of 
the squares error of 15.26. The values obtained for the 
variables were bu=2.7883, b225=0.8934, bjjsO.OOOl and b42=0.1778. 
Again the value for as<, dropped to the lower limit allowed by 
the constrained optimizer, 0.0001, and remained there. Figure 
3.10 illustrates the time-temperature profiler of the node 
three match. 
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T-gure 3.10 Time-Teinperature Profiles Node Three Match 

(3) Node Four Temperature Modeling, Again the 
program was modified to match the experimental data for node 
four only and the simulation converged for this set up with a 
sum of the squares error of 39.539. The values obtained for 
the unknowns were bus2.92, bjjal.SS, bj,*0.212 and b 4 i= 0.02 and 
again the value of a^c went to the lower boundary limit sec in 
the constrained optimizer and remained there. The time- 
temperature profiles are illustrated in Figure 3.11. 
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Figure 3.11 Time-Tenperature Profiles Node Four Match 

(4) Results. In all three simulations, using only 
a single set of experimental data, the program converged and 
the calculated and experimental profiles of the node that was 
being matched were as accurate as the optmizer could get. The 
error using node two data was less than that of node three and 
much less than that of node four. The problem arises in the 
calculated temperature histories of the other two nodes not 
being matched. It can be seen in Figures 3.9, 3.10 and 3.11 
that the two nodes not being matched are significantly 
different from the experimental data. In the simulation using 
the experimental data from node two the profiles of node two 
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and three match fairly closely while the profile for node four 
is nowhere near the experimental data. In the simulation 
using the node three experimental data, node three profiles 
have a fair match but the profiles for nodes two and four are 
unsatisfactory. The same condition holds when using the 
experimental data for node four only. It must be noted that 
the magnitude of the error between the calculated and 
experimental data is dependent on the quality of the 
experimental data. If the data follows a first order linear 
response closely then the match of the calculated and 
experimental profiles is adequate and the sum of the squares 
error is low. If the experimental profiles vary from an 
exponential curve then the first order linear model matches 
the experimental profile the best it can and the error will be 
large. 

o. Huiuml Adjustmmt of Vmriablma-Fomrd Modal 

The results of the single node tenperature profile 
simulation permit the use of the forward model to generate a 
series of time-temperature profiles with the final values 
obtained for the unknowns. By adjusting the various unknowns 
we can see the effect one variable has on another. A set of 
values needs to be selected for the unknowns. These values 
can be arbitrary since we are looking at the effects of one 
value on another. The initial starting values for the 
unknowns were chosen as, buai.ve, b2as0.0455, bj 2 => 0 . 06 , b 4 i= 0.02 
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and aso-O.OOOl. The value of aso was chosen as 0.0001 because 
in all the previous simulations this unknown juitped to this 
lower value limit and remained there. The values for the 
other unknowns were chosen by scaling up the final values of 
the quarter scale model results obtained by Reno [Ref. 1] and 
considering the values obtained from the single node 
teirperature profile results. The effects of the variables 
were that, the value of bu influenced the teirperature profile 
of node two and had little or no effect on nodes three and 
four, while the same holds true for the other convective rate 
coefficients. Altering each one consecutively varied the 
result of that particular node profile but had little effect 
on the outcome of the other node calculated temperature 
profiles. As an example of this effect the values of the 
unknowns listed above were used in the forward model. The 
model used a set of recovery tenperatures with values of 
Tri=2970 and Trj=2870. There was no need in normalizing these 
recovery tenperatures in the forward model. The effects of 
one variable on the other was being examined and not actually 
trying to determine the proper values for the unknowns as done 
in the PID modeling program. The profiles generated by this 
forward model simulation are illustrated in Figure 3.12 where 
the three nodes of concern, nodes two, three, and four are 
displayed. 
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Figure 3.12 Forward Model Time-Tcnperature Profiles 
These profiles show a good match between the calculated and 
experimental tenqperature profiles of node two. The profiles 
of nodes three and four do not follow the experimental data 
but again the experimental data for these two nodes does not 
follow an exponential curve. Another thing that needs to be 
pointed out is the calculated curves for the three nodes 
appear to be linear straight lines. This is not an 
exponential curve as esqpected from a first order linear model. 
To show the effects of changing one variable on the other 
temperature profiles, the value of b)} was changed to 0.03 and 
the forward model simulation was executed. The resulting 
time^temperature profiles are displayed in Figure 3.13. 
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Figure 3.13 Forward Xodel Slaulatlon Node 3 Variant 
Comparison o£ Figures 3.12 and 3.13 show that there was a 
significant change in the calculated tenperature profile of 
node three, while the profiles of nodes two and four changed 
very little in magnitude cou^red to node three. The program 
calculated the sum of the squares error for both simulations. 
The error for the first simulation was 74.4 while the error 
for the second simulation was 213.577. 

Adjustments can be made to each of the unknown 
variables in the forward model to create a better match for 
any of the individual curves. This process of adjusting the 
unknowns is a manual atten^t to do what the optimizer was 
intended to do in the PIO program. The manual adjustment 
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allows the user to view the effects that individual rate 
coefficients have on the temperature profiles of the other 
nodes. 

d. Ttio Node T&uperatiire Profile Modeling 

The previous section established the convergence of 
the modeling program and the matching of a single time- 
temperature profile. The model will be upgraded to matching 
two temperature profiles in an atten^t to obtain a three node 
full scale matching model. 

(1) Nodes Two and Three Temperature Modeling. The 
modeling program was modified to match the temperature 
profiles of two nodes, vice one, with the sum of the squares 
error being calculated based on both temperature profiles. 
The first simulation will attenpt to match the temperature 
profiles of nodes two and three. The assumptions for the two 
node matching simulations are the same as those for the single 
node in that all the unknowns are considered independent of 
each other and the constrained optimizer will be utilized with 
the upper boundary of 100 and a lower boundary of 0.0001. The 
simulation converged with the resulting sum of the squares 
error of 15.217 and values for the unknowns bii-2.83, 
b22**0.0297, b32'"0.0056, b42'<‘0.9408 and asg-O.OOOl. Again the 
value for asg jumped to the lower limit and remained there 
while the optimizer was changing the values for the other 
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unknowns. The time-temperature profiles for the two node 
matching process are illustrated in Figure 3.14. 



Figure 3*14 Two Hod* Tcoporaturo Profllos 
The simulation program matched the time-ten^erature profiles 
for nodes two and three while the calculated profile for node 
four was no where near the experimental profile. Figure 3.14 
illustrates that the calculated profile for node four, T4 
calculated profile, was unsatisfactory. 

(2) Nodes Two and Four Tenperature Modeling. The 
program was again modified to match the data for the 
temperature profiles from nodes two and four. The simulation 
using these experimental profiles yielded a sum of the squares 
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error of 216.02 and values of the unknowns 1.2010, 
b22-0.1151, b32-0.0001, b42-0.0001 and agQ-O.OOOl. The time 
temperature profiles for this simulation are illustrated in 
Figure 3.15. 



Flgiure 3.15 Ti«o Node Taosperature Prof lias 
The experimental and calculated profiles for node two matched 
adequately, but the calculated profile for node four is 
substantially off from the experimental data. The value of 
® 5 G reached the lower limit set by the optimizer and 

remained there. 

(3) Nodes Three and Four Temperature Modeling. The 
final two node matching simulation was conducted by modifying 
the main program and calculating the sum of the squares error 
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based on the temperature profiles from nodes three and four. 
The simulation yielded values of bii=2.19, b22«0.77, b32»0.06, 
b42»0.0001 and agg-O.OOOl with a sum of the squares error of 
185.95. As can be seen in Figure 3.16 the program adequately 
matched the profiles for nodes three and four but there was a 
significant difference in the calculated and experimental time 
temperature profile of node two. 



(4) j?esults. The program was able to match the 
ten^erature profiles from node two and to a lesser extent node 
three, but had difficulty in matching the profiles for node 
four. This is due to the fact that the experimental profile 
for node four is not an exponential curve. The first order 
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linear modeling program provides an descent attempt to model 
the experimental data and provides as close an estimation 
possible. The use of the node four experimental data causes 
a significant increase in the sum of the squares error. 
Although it may be possible to get an adequate match for two 
temperature profiles the third is by far not satisfactory, 
e. Convergent Simulation Model 

(1) Three Node Temperature Profile Matching. It 
has been demonstrated that the modeling program can match one 
and two temperature profiles. The process is now upgraded to 
model all three ten^erature profiles of the vane, that is the 
modeling of nodes two, three, and four. This is the primary 
reason behind the modeling of the full scale vane. 

(2) Unknowns and ABauaptions. The model simulation 

for the three node ten^erature matching will assume that the 
unknown variables, b^]^, b 22 « ^321 b 42 and a^Q, are independent 
of each other. Previous five node model simulations 
considered a constant in the same manner as the quarter 
scale modeling did. Changing the value of in the quarter 

scale model had little effect on the values of the other 
imknowns . This fact has not been established in the full 
scale modeling process which was the reasoning behind nvaking 
bj^]^ a variable. The initial conditions will be set to zero 
for all five unknowns and the values for the recovery 
temperatures are T}y^<*i2670 and Tr 2*2570 Kelvin. For this 
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simulation the sxim of the squares error were be calculated 
based on all three temperature profiles and the boundary 
limits for the constrained optimizer were set to 0.0001 for 
the lower boundary and 100 for the upper boundary. 

(3) Results. The simulation was executed for a 
time period of zero to five seconds. The simulation converged 
with a sum of the squares error of 58.5 and the values 
obtained for the unknown variables were bj^j^»1.3787, b22»0.0938, 
b32"0.0862, b42»0.0484 and asQ-O.OOOl. As with the single and 
two node temperature matching the value for agQ dropped to the 
lower boundary set in the optimizer. The low value obtained 
for a 5 g relates to the negative value obtained for the similar 
unknown a 4 Q of the quarter scale model from Oriels [Ref. 3] 
simulation. The constrained optimizer prevented the value for 
a^g frc»n going negative in the five node model simulation. 
The tlme'ten^erature profiles generated by this simulation are 
illustrated in Figure 3.17. 




TIME (sec) 


Ftgur« 3.17 fiv* Hod* Van* ICodsl Taap^rature Blstorlta 

The match between the profile labeled T2, for node two, was an 
adequate match as it was in some of the previous simulations. 
The match between the experimental and calculated profiles for 
nodes three and four were not as close as node two. Again 
this is due to the quality of the esqterimental data taken. If 
the data does not follow am exponential rise such as a first 
order linear system, then the optimizer is less accurate in 
matching the experimental temperature profiles. 

The simulation above was conducted using the recovery 
temperatures of Tri« 2670 and Tg2-2570 Kelvin. These values 
were obtained from past work conducted by Nunn (Ref. l] and 
Reno (Ref. 2J. Since there was no direct way to determine 
these values for the full scale oxperiiuental rocket motor 
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firings, due to the unavailability of firing data, the values 
of the recovery temperatures were altered to illustrate the 
effects on the model simulation. The first modification was 
to change the values of the recovery temperatures to the value 
of the chamber temperature. This would make Tj^j^=2970 and 
Ti^ 2“2870 degrees Kelvin. This simulation converged with a sum 
of the squares error of 56.5, The values of the unknowns were 
bii=1.4747, b22=0.0716, b32=0.0761, b42=0.0430 and a5G=0.0001. 
The time-temperature profiles were very similar to those in 
Figure 3.17. Another simulation was conducted with the values 
of Tjjj^=3550 and Tjj2“3450 degrees Kelvin. This simulation 
yielded an error of 54.5 and values of bj^j^=*2.3045, b22"0*033, 
b32“0.0621, b42«0.0354 and aggaO.OOOl. The time-temperature 
profiles for this simulation are again very similar to those 
obtained in Figure 3.17. This would be expected due to such 
a small variation in the sum of the squares error. 

f. Simulation With various Initial Conditions 

To test the robustness of the simulation program 
the initial conditions were varied. Changing the initial 
conditions changes the starting point for the model 
simulation. Variation of the starting values for the unknowns 
yielded the same final values. The program converged to the 
same point and the same profiles when the initial conditions 
were changed to different values. 
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g. Calculations of Heat Transfer Coefficients 

The values obtained for the unknovm variables, for 
the three node temperature matching simulation with the 
recovery tenperatures of Tj^3^=2670 and Tj^=2570 Kelvin, were 
b2^^®1.3787, b22’*0*052®i b 22 “ 0 • , b^ 2 *®• and agQ=0.0001« 

The convecuive heat transfer coefficients can now be 
determined but the values for the resistance must first be 
determined. These calculations yield Rjf3^=0.1044 K/W, 
Rp22-0.3101 K/W, Rp23”0.2530 K/W and Rp24-0.3605 K/W with the 
values for the corresponding convective heat transfer 
coefficients; ht-22027.5 W/ro^*K, W/m^-K, h^3=l057 

W/m^*K, and hv4"594 W/m^«K. Where hj. is the convective heat 
transfer coefficient for the tip and 1 x^ 2 » ^^3 Nr 4 
convective heat transfer coefficients from the free stream to 
nodes two, three and four of the vane, respectively. 

The simulation involving the use of Tr 3^«2970 and Tr 2“2870 
Kelvin ^/ielded the values of bj^;^-1.4747, b22“0,07l6, 

b32-d.076l, b42«0.0430 and a^Q-O.OOOl. These values yield 
resistances of Rpi-0.0976 K/W, Rp22-0>4062 K/W, Rp23-0.287 K/W, 
and Rp24-0.4058 K/W with the resulting convective heat 
transfer coefficients of h^«2356l.3 W/m**K, w/m^*K, 

1x^3-932.86 W/m^*K and hv4-527.74 W/m^*K. The last simulation 
involved changing the recovery temperatures to "^Rl -3550 and 
Tr 2 -34*30 Kelvin. The results of this simulation yielded 
b2i,-2, 1045, b22"0*033, b32*0.C621, b42»0.03S4 and asg-0.0001. 
The rc««i8tance values obtained frcoa these results are 
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Rfi=0.0624 K/W, Rf22=0.8815 K/W, Rf 23=0.3512 K/W and Rf 24=0.4928 
K/W. Using these resistances the values for the convective 
heat transfer coefficients are ht«36819 W/m^^K, hy 2 = 404-6 
W/m2*K, h^3=761.5 W/m2«K and h^4=434.5 W/m^^K. 
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DISCUSSION OF RESULTS 


The full scale six noae model atteir^Jted to directly scale 
up the four node quarter scale vane model. The decision to 
add the two extra nodes was necessary in order to model the 
larger vane and was believed to more accurately describe the 
heat transfer process for the larger amount of material 
contained in the full scale vane. The first series of 
simulations involved the use of the 2XSSQ optmizer to adjust 
the unknowns in order to obtain a match between the 
experimental and calculated time-temperature profiles. The 
unknown rate coefficients for the six node model were 
identified as b^,, b} 3 « b^, as(, a^s and aso. The assumptions 
were that b^ was a constant, as it was in the quarter scale 
simulations, and that the convective rate coefficients, bjj, 
bj] and bt^, were equal. These series of simulations failed to 
converge due to the substantial variation in the unknown 
values, Ssc a^s and a«, which contain the resistances between 
the shaft and mount and the inount and ground. The values of 
a^s and a^a also contain the unknown value of the mount node 
capacitance. During the simulation these variables jumped 
substantially in value causing the differential equations to 
become stiff. The juiqp in the variables were both in the 
negative and positive direction. The experimental data for 
these simulations was taken from the shaft and mount locations 
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in the same manner that the experimental profiles were Ui;;ed in 
the quarter scale modeling process. When the simulation 
failed to converge the experimental data was changed to that 
relating to nodes two, three and four of the six node model. 
In the physical test firings, for the full scale vane, the 
tenperature profiles corresponding to the shaft and mount 
locations had minimal variation from the ambient. It was 
believed that the simulation program had difficulty adjusting 
to the small changes in the mount and shaft temperature 
profiles while trying to make large adjustments to the 
unknowns from nodes two, three and four. That is, over a time 
increment the tenperature change was small in the shaft and 
mount but was large in the tip and fin area. This idea led to 
abandoning the uf' of the experimental data from the shaft and 
mount and using the experimental data from the fin area nodes, 
nodes two, tnree and four. The simulation using this new set 
of experimental profiles also failed to converge with the same 
results for the unknowns, a^^ and a(o> 

^vjther idea considered controlling the simulation by 
preventing the unknown values from varying substantially. 
This was achieved by the use of the constrained optimizer, 
which would permit the setting of boundary limits to hold the 
unknowns within a certain range. This range was selected as 
positive values with a lower boundary limit of 0.0001 and an 
upper boundary limit of 100. The assumptions for the unknowns 
remained the same as the previous simulations. This 
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simulation using the constrained optimizer also failed to 
converge, where the values of the un^novms asg, a^s and agg 
reached the set limits and remained there. Attenpts were made 
at raising the upper boundary in the hope that the higher 
boundary limit would allow the optimizer a chance to avoid or 
recover from any local minimum and continue on with the 
simulation. Raising the upper boundary limit had the same 
result as before, the unknown values reached the boundary 
limits and remained there and the program failed to converge. 
Any attenpt at setting the lower bou.ndary limit to a negative 
value or raising the upper limit to too high a value caused an 
arithmetic fault to occur. A set of simulations was conducted 
using each optimizer but changing the assumptions, making the 
convective rate coefficients bjj, bj), and b^j independent of 
each other. In both cases these simulations failed to 
converge. 

Failure to resolve the problem of the unknowns associated 
with the shaft and mount led to the modification of the model 
from six nodes to five. This was done by eliminating the 
mount node which eliminated the unknowns as«, a«s and a^o and 
added the unknown ajo. This was justified by the fact that 
there was no consistent mount arrangement used from one test 
firing to another. The unknown heat transfer coefficients 
that are being determined are from the free stream to the fin 
area of the vane, not the mount or shaft. The small 
temperature change in the experimental profiles indicated that 
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very little energy from the fin section was being transmitted 
to the shaft and mount. The unknowns were identified as bn, 
b22/ b32, b42 and The assun^Jtions, for the first series of 
simulations, were that b^ was a constant and b22=b32=b42. With 
the program modification in place the simulation was conducted 
using the ZXSSQ optimizer. The result of this simulation was 
the same as in the six node model, the program failed to 
converge using this optimizer. Attenpts were made at using 
only a single experimental temperature profile, but these 
simulations also failed to converge. The five node program 
was modified as was the six node to use a constrained 
optimizer to place a set of boundary limits on the unknowns 
with the limits set to the same values as in the six node 
simulations. The assumptions for the five node model unknowns 
remained the same with constant and b22=b32=b43. This 
simulation yielded the same results as the previous modeling 
attempts in that it failed to converge. A more systematic 
approach was taken to determine if the simulations could be 
made to converge. To accomplish this the simulations would be 
constructed to match only one temperature profile at a time. 
The assumptions were that none of the unknowns were dependent 
on one another and the term bj, was no longer considered a 
constant as before. The unknown variables were identified as 
^11 ♦ ^42 and a^g. 

The simulations converged and succeeded in providing an 
adequate match for the selected single node time ten^erature 
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profiles. The other two fin area node profiles were 
unsatisfactory in the match between the experimental and 
calculated profiles. The model was modified to match two 
tenperature profiles from the fin area of the vane. This 
series of simulations converged and an adequate match was 
achieved for two of the three tenperature profiles, but again 
the problem was in the third profile. The program matched the 
two temperature profiles that it was set up to match but the 
calculated profile for the third node was substantially off 
from the experimental data. Being able to match two node 
temperature profiles the program was modified in an attempt to 
match all three nodal temperature profiles. The program was 
modified again to read experimental data for nodes two, three 
and four, and to match all three time temperature profiles. 
This simulation was conducted using the assumptions that b^, 
baa# baa# b4a# and aso were independent, the recovery 
temperatures were Tni =2670 and =2570 Kelvin, and the initial 
conditions set to zero. This simulation converged and yielded 
the profiles displayed in Figure 3 .. 17 . The simulation was 
tested for its robustness by starting with various initial 
conditions. The program converged to the same set of values 
with the same set of time-temperature profiles. The match 
between the calculated and experimental profiles of node two 
were satisfactory while the match for nodes three and four 
could only be considered adequate. The deficiency here is not 
with the modeling program but with the experimental time- 
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teirperature profiles. The experimental data does not follow 
an exponential curve, where as the modeling program is a first 
order linear model and can only match exactly an exponential 
profile. The modeling program matches these profiles the 
closest it can to an exponential curve. Another point to note 
is the calculated teirperature profiles for all three nodes are 
almost a straight line with very little curve in the 
calculated values in any of these three temperature profiles. 
Since the exact recovery tenperatures were unknown the three 
node temperature matching simulation was conducted using two 
other sets of profiles. Table 4 lists the recovery 
tenperatures used and the corresponding convective heat 
transfer coefficients. 

TABLE 4 


Recovery 

Temperatures 

T.,/T„ 

Convective Heat Transfer Coefficient 
W/m*»K 

ht h„3 Ki Ki 

2670/2570 

2970/2870 

3550/3450 

22027 1150 1057 594 
23562 878 933 527 
36819 405 762 434 


The values of the heat transfer coefficients displayed in 
Table 4 are of the same magnitude as those obtained by Oriels 
in the quarter scale simulations. The values for the heat 
transfer coefficients from quarter scale modeling were 
heB30122 W/m^*K and hy=210 W/m*«K (Ref. 3:p. 12]. 
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The calculated profiles in the convergent five node model 
simulation were nearly straight lines. This effect was caused 
by the reduction of the conductive rate coefficients, the 
values. These values are determined from equations (2) and 
(4), the resistance and capacitance equations. Consider the 
area A as a length squared and the volume V as a length cubed. 
The value of is equal to (RC)*S and inserting the 
equivalent values from equations (2) and (4), and considering 
the physical properties as constants, the value of RC becomes 
a function of the length squared. This means that the value 
of ajj is a function of the (length^)■^. As the quarter scale 
model is increased in size the conductive coefficients aij 
decrease as the reciprocal of the length squared, if the same 
idea is applied to the convective coefficients it will be seen 
that the value of RC changes as the length changes and the 
convective rate coefficient, bij, changes as the (length) 

As the model goes to full scale, the geometric properties 
become larger and thus the rate coefficients loecome smaller. 
Since the conductive values decrease as the reciprocal of the 
square of a length, where the convective values only decrease 
as the reciprocal of the length, the convective rate 
coefficients become dominant. 

The quarter scale modeling conducted by Oriels (Ref. 3], 
demonstrated that the value of the convective rate 
coefficient, b^, from the tip section of the vane could be 
adjusted with minimal effect on the other unknown variables. 
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This was due to the ablative cooling of the tip due to the 
erosion of the tip area of the vane by the rocket exhaust. In 
the full scale experimental rocket firings there was minimal 
vane erosion as cortqpared to the quarter scale vane tests. The 
energy in the full scale tests was not being carried away by 
the ablative cooling but was being transmitted into the vane. 
This was one of the primary reasons that led to making the 
value of the convective rate coefficient, bn, a variable in 
the full scale model instead of the constant that was assumed 
based on the quarter scale model results. 
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CONCLUSIONS 


• The four node quarter scale model cannot be directly 
scaled up to a full scale model due to the problems 
encountered with substantial juii¥>s in the unknowns asg and 
ags, corresponding to the shaft and mount locations. 

• The use of a constrained optimizer is necessary to hold 
the value of the unknowns to positive values and prevent 
the equation solver from experiencing an arithmetic fault 
or causing a fatal error because the energy balance 
equations became stiff. 

• The value obtained in the quarter scale vane modeling 
process more adequately reflects the average convective 
heat transfer coefficient of the fin area where as the 
v,^l*aes obtained from the full scale modeling display a 
3iore accurate representation of the heat transfer process 
in the vane. 

• The calculated profiles of Figure 3.17 indicate that as 
the vane size increases the convective heat transfer 
properties become the dominant player in the heat transfer 
process of the vane. 

• The rise in the shaft and mount experimental temperature 
profiles is likely due to the rocket exhaust plume 
impinging directly on the shaft of the vane. 

• The value of the convective rate coefficient, bj,, must be 
considered an unknovm variable in the time*temperature 
matching process, due to the minimal amount of vane 
erosion and the lessening effect of ablative cooling in 
the full scale vane. 
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Al^C^S VCR FORIHER STUDY 


REC< 


• It is recommended that further study modify the five node 
model by reducing it to three nodes and coir^aring the 
results to see if a single node in the fin area of the 
vane could produce an acceptable convective heat transfer 
coefficient. 

• The full scale model needs to be able to take into account 
the effects of radiation from the casing of the rocksst 
firing chamber. Itiis may condensate for the non- 
linearities in the experimental temperature profiles. 

• Vane erosion was evident in the quarter scale testing and 
some vane erosion occurred in the full scale testing. The 
model needs to be upgraded to include a conponent which 
can compensate for the erosion based on the metal content 
of the rocket exhaust. 

• Since the convective rate coefficient cannot be considered 
a constant at the tip, a thermocouple attached to the tip 
of the vane to collect the experimental time-tenperature 
data would greatly enliance the thermal model develcjcanent. 

• The process for the attachment of theosKJCouples needs to 
be addressed. The experimental data had glitches or 
spikes in the temperature profiles during various times of 
the rocket firing indicating that the therawcouple 
temporarily lost contact with the vana. 
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APPEHDIX A. PIO PROGRAM 

The progrcun contained is the FORUL^N code main program 
used in the vane modeling process. Hie programs name was 
changed to NODES to distinguish it from the PID program using 
the six node modeling approach 
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Ri *r> Ri If Jf> 


Program NODES 

This program is the PID program for the five node vane model, 
external temp 

integer m,n,iparm(6),ibtype,ldfjac 
parameter (m=183,n=5,ldfjac=m) 
real*8 rparm(7),x(n),f(ra),xjac(m,n),xg(n), 
xlb(n),xub(n),xscale(n),fscale(m), 
ssq,float, 

a(5,5),b(5,5),u(5),t2(61),t3(61),t4(61),ys(5,61), 
rho,k,op,vt,vf,vs,atf,af s,Itf,If s,sf, 
atf1,atf2,atf3,vf1,vf2,vf3,ubl,ub2 
intrinsic float 

common/datal/a,b,u,t2,t3,t4,ys 

Open files for data input/output 

open(10,name= • result.dat •, status= • new') 
open(9,name= ' temp.mat ', status='new') 
open(8,name='datarn.dat ‘, status® ^ old') 
open(7,name® • input.dat', status®'old') 

read in experimental data 

do i®l,61 

read(8,*) t2(i) 

enddo 
do i®l,61 

read(8,*) t3(i) 

enddo 
do i®l,61 

read(8,*) t4(i) 

enddo 
close(8) 

read in input data 


rsad(7,*) 
read(7,*) 
rcad(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
close(7) 


rho,k,cp 


vt, vfl, vf2, vf3, vs 


atfl, atf2, atf3, afs 

Itf, ifs 


sf, ubl, ub2 
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full scale data 


rl2=^100. o*ltf/ (k*atf 1) 
r23=100.0*ltf/(k*atf2) 
r34=100.0*ltf/(k*atf3) 
r45=100.0*lfs/(k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*cp*vf1*0.000001 
c3=rho*cp*vf2*0.000001 
c4=rho*cp*vf3*0.000001 
c5=rho*cp*vs*0.000001 


scaled data 

rl2=rl2/sf 

r23=r23/sf 

r34=r34/sf 

r45=r45/sf 

cl«cl*sf**3 

c2=c2*sf**3 

C3=c3*sf**3 

C4=c4*sf**3 

c5=c5*sf**3 

rfl=0.5 

do 

u(i)=0.0 

do j*l/5 
a(i,j)»o.o 
b(i, j)=*0.0 
enddo 

enddo 

a(l,2)-l/(cl*rl2) 

a(2,l)-l/(c2*rl2) 

a(2,3)»l/(c2*r23) 

a(3,2)-l/(c3*r23) 

a(3,4)"l/(c3*r34) 

a(4,3)“l/(c4*r34) 

a(4,5)»l/(c4*r45) 

a(5,4)»l/(c5*r45) 

b(l,l)“l/(cl*rfl) 

wirite(6,*)^ al2 a21 a23 a32' 

writet6,3002) a(l,2),a(2,1),a(2,3),a(3,2) 
write(6,*)' a34 a43 a45 a54' 

write(6,9003) a(3,4),a(4,3),a(4,5),a(5,4) 

002 format(5f9.3) 

303 format(4f9.3) 

a5g "O.o 
b(l,i)-o.o 

b(2|2)«0.0 
b(3,2)«0.0 
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a(l,l)=-(a(l,2)+b(l,l)) 
a(2,2)=-(a(2,l)+a(2,3)+b(2,2)) 
a(3,3)=-(a(3,2)+a(3,4)+b(3,2)) 
a(4,4)=-(a(4,3)+a(4,5)+b(4,2)) 
a(5,5)=-(a(5,4)+a5g) 

u(l)=ubl 

u(2)=ub2 

xg(l)=a5g 

xg(2)=b(l,l) 

xg(3)=b(2,2) 

xg(4)=b(3,2) 

xg(5)=b(4,2) 

set up parameters for DBCLSF call 
do lc=l,n 
xscale(]c)=1.0 
xlb(k)«0.0001 
end do 

xub(l)=100.0 
xub(2)=100.0 
xub(3)=100.0 
xub(4)=100.0 
xub(5)=100.0 

do 1=1,m 
fscale(l)»l.0 
end do 

ibtype»0 


call dbclsf(temp,m,n,xg,ibtype,xlb,xub,xscale , fscale , 
& iparm,rparm,x,f,xjac,ldfjac) 

calculate unknown resistances and capacitances 

aSg *x(l) 
b(l,l)-x{2) 
b(2,2)«x(3) 
b(3,2)-x(4) 
b(4,2)=x(5) 


rfl •l/(b(l,l)*cl) 
rf22-l/(b(2,2)*c2) 
rf23-l/(b(3,2)*c3) 
rf24»l/(b(4,2)*^C4) 
r5g«l/(a5g*c5) 

print and save results 

write(6,*) • a5g bll b22 b32' 

write(6,9009) x(l),x(2),x(3),x(4),x(5) 


09 format(5fl2.2) 
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;UU 


rormat;(5ti2.4) 


write(l0,*)'a5g b(l,l) b(2,2) b(3,2)' 

write(10,9000) x(l),x(2),x(3),x(4),x(f) 
write(10,*) 
write(10,*) 

write(10,*)'rfl rf22 rf23 rf24 r5g' 

write(10,9000) rfIrf22,rf23,rf24,r5g 

write the temp-time data for MATLAB analysis 

do i=l,61 
tt=0.0768*float(i) 

write(9,9001)tt,ys(2,i) ,ys(3,i) ,ys(4,i),t2(i),t3(i),t4(i) 
enddo 

01 format(lf6.2,6fl0.3) 

close(lO) 

close(9) 

end 


Subroutine TEMP (m,n,x,f) 

This calculates the temperature-time history using the current 
parameters supplied by ZXSSQ called from PID. It calculates 
an error function returned to ZXSSQ based on the differences 
between predicted and observed temperature histories. 

integer maxparam, neq 
parameter(maxparam^SO, neq=5) 
integer idO#istep,nout,m,n 

real*8 t,y(neq),a(neq,neq),b(neq,neq),u(neq),param(raaxparam) 
real*8 tend,tol,fen,float,a5g,ys(5,61),x(n), 

& f(m),t2(61),t3(61),t4(61) 

intrinsic float 
external fcn,divprk,sset 

common/datal/a,b,u,t2,t3,t4,ys 

open (12, name® • incoming. dat •, .status®' new •) 

aSg »x(l) 
b(l,l)»x(2) 
b(2,2)-x(3) 
b(3,2)«x(4) 
b(4,2)»x(5) 

write(6,8000)a5g,b(l,l),b(2,2),b(3,2),b(4,2) 

00 format(5fl2.4) 

a(l,l) —(a(l,2)+b(l,l)) 
a\'2,2) —(a(2,l)+a(2,3)+b(2,2)) 
a(3,3) —(a(3,2)+a(3,4)+b(3,2)) 
a(4,4) —(a(4,3)+a(4,5)+b(4,2)) 
a(5,5) —(a(5,4)+a5g) 
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t=0.0 

do i=l,neq 
y(i)=»0.0 

do j=l,6l 
ys(i,j)=0.0 
enddo 

enddo 

tol=0.0005 

call sset (maxparm, 0.0« param, 1} 


idO=l 

do isteps'ljei 

tend=0.0768*float(istep) 

CALL DIVPRK (idO, neq, fen, t, tend, tol, param, y) 

do i=l,5 

ys(i,istep)=y(i) 

enddo 

enddo 

Final call to release workspace 
id0=3 

call divprk (idO,neq,fen,t,tend,tol,parain,y) 

calculate error functions 

do i=l,61 
f(i)-ys(2,i)-t2(i) 
f(i+61)-y8(3,i)-t3(i.) 
f(i+l22)»y8(4,i)-t4(i) 
enddo 

print out rms error 
ssqr»o.0 
do i-1,122 
ssqr»ssqr+f(i)* f(i) 
enddo 

ssqr»ssqr/122 

xer“sqrt(ssqr) 

write(6,*) xer 

write(12,*) xer 

return 

end 


subroutine fcn(neq,t,y,yprime) 
integer neq 

real*8 t,y(neq),yprime(neg) 
real*8 a(5,5),b(5,5),u(5),d,ys(5,61) 
conmon/datal/a,b,u,t2,t3,t4,ys 

thrust profile simulation as step input 


if (t.gt.0.2) then 
d*1.0 
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else 

d=0.0 
end if 

do i=l,neq 
ypriine(i) =0.0 

do j=l,neq 

ypriine(i)=ypriine(i)+a(i, j)*y(j)+b(i, j)*u(j) *d 
enddo 

enddo 

return 

end 
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APPENDIX B. GEOMETRIC AND PHYSICAL DATA FILE 


The file in this appendix, INPUT.DAT, is the 
containing the geometric and physical properties of the 
The values are listed in full scale. 


file 

vane. 
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NAWC INVERSE HEAT TRANSFER PROGRAM. INPUT DATA FOR FULL SCALE 

Material properties: rho (kg/in''3), k (w/m-deg k), Cp (J/kg deg k) 

18310.0, 173.0, 146.0 

Vol (tip, cin*3), Vol (fin, cm*3), Vol (shaft, cm^3) 

2.6 12.86 17.15 21.44 23.0 

X-section areas: tip-fin (cm''2) fin-shaft (cm''2) 

5.6 5.9 6.2 5.2 

Conductive path tip-fin (cm) fin-shaft (cm) 

2.774 5.0 

Scale factor: 

1.00 2670 2570 
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APPENDIX C. TPROFILE - FORNARD MODELING PROGRAM 
The program this appendix is the forward modeling program 
based on a modification of the main program, PID, and used to 
illustrate the effects of the unknown variables on the 
resultant time tenperature profiles. 
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ta 


Program TPROFILE 

This is a forward model program developed from the program 
PID by removing the optimizer. 

external temp 

integer m,n,ixjac,nsig,maxfn,iopt,infer,ier 
parameter (m=183,n^S) 

real*8 parm(6),x(n),f(m),xjac(m,n),xjtj{(n+1)*n/2), 
work(5*n+2*m+((n+l)*n/2 )), 
eps,de1ta,ssg , ubl,ub2, 

a(5,5),b(5,5) ,U(5) ,t2(61),t3(61),t4(61),ys(5,6?)/ 
rho,k,cp,vt,vf,vs,atf,afs,Itf,lfs,sf 
common/datal/a,b,u,t2,t3,t4,ys 

Open files for data input/output 


open(10,name-^result.dat *, status='new') 
open (9, names* • temp. mat •, status* • new •) 
open(8,name* * datam.dat', status*'old') 
open(7,name* • input.dat•, status* • old') 

do i*l,6l 

read(8,*) t2(i) 

enddo 
do i“l,61 

read(8,*) t3(i) 

enddo 
do i*l,61 

read(8,*) t4(i) 

enddo 
close(8) 

read in input data 


read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read (7, 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
read(7,*) 
cloae(7) 


rho,k,op 

vt, vfl, vf2, vf3, vs 
atfl, atf2, atf3, afs 
Itf, Ifs 
sf, ubl, ub2 


initial conditions 


Cull scale data 
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rl2=l00.0*ltf/(k*atf1) 
r23=100.0*ltf/(k*atf2) 
r34=l00.o*ltf/(k*atf3) 
r45=100.0*lfs/(k*afs) 
cl=rho*cp*vt*0.000001 
c2=rho*cp*vf1*0.000001 
c3=rho*cp*vf2*0.000001 
c4=rho*cp*vf3*0.000001 
c5=rho*cp*vs*0.000001 

write(6 ^ *)r12,r23,r34,r45,cl,c2 ^ c3,c4,c5 

scaled data 

rl2=rl2/sf 

r23=r23/sf 

r34=r34/sf 

r45=r45/sf 

cl=cl*sf**3 

c2=c2*sf**3 

c3=c3*sf**3 

c4®c4*sf**3 

C5»c5*sf**3 

rf1-0.076318 

do i-1,5 
u(i)»0.0 

do j-l,5 
a(i,j)=0.0 
b(i,j)»0.0 
enddo 

anddo 

a(l,2)-l/(cl*rl2) 

a(2,l)-l/(c2*rl2) 

a{2,3)-l/(c2*r23) 

a{3,2)»l/(c3*r23) 

a(3,4)-l/(c3*r34) 

a(4,3)-l/(C4*r34) 

a(4,5)-l/(c4*r45) 

a(5,4)«l/(c5*r45) 

b(l,l)»j./(cl*rfl) 

write(6,*)' al2 a2l a23 a32 bll' 

write(6,9002) a(l,2),a(2,l),a(2,3),a(3,2),b(l,l) 
write(6,*)' a34 a43 a45 a54' 

write(6,9003) a(3,4),a(4.3),a(4,5),a(5,4) 

)02 format(5f9.3) 

)03 Corinat(4f9.3) 

aSg -0.0001 

b(l,l)-300.0 

b(2,2)-1.339 

b(3,2)-0.56S7 

b(4,2)-0.0037 
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a\j., J./ -vcna.,6/-rw^x,x;; 

a(2,2)=-(a(2,l)+a(2,3)+b(2,2) ) 
a(3,3)=-(a(3,2)+a(3,4)+b(3,2)) 
a(4,4)=-(a(4,3)+a(4,5)+b(4,2)) 
a(5,5)=-(a(5,4)+a5g) 

u(l)=ubl 

u(2)=ub2 

x(l)=a5g 

x(2)»b(l,l) 

x(3)=b(2,2) 

x(4)=b(3,2) 

x(5)=b(4,2) 

call teinp(x,m,n,f) 


wrice the temp>time data for MATLAB analysis 
do 

tt».0768*float(i) 

write{9,900l)tt,ys(2,i),ys(3,i),ys(4,i),t2(i),t3(i),t4(i) 
enddo 

01 format(If6.2»6f10.3) 

close(10) 
close(9) 
end 


Subroutine TEMP 

This calculates the temperature-time history using the current 
parameters supplied by ZXSSQ called from PID. It calculates 
an error function returned to ZXSSQ based on the differences 
between predicted and observed temperature histories. 

integer maxparam, neq 
parameter(maxparam^SO, neq«S) 
integer idO,istep«nout^m,n 

real*8 t,y(neq),a(neq,neq}«b(neq,neq),u(neq),param(maxparam) 
real*8 tend,toi,fcn,float,a5g,ys{ 5,61} ,x(n), 

& f(m),t2(61),t3(61),t4(61) 

intrinsic float 
external fen,divprk,sset 

common/datal/a,b,u,t2,t3,t4,ys 

aSg -x(l) 
b(l,l)-x(2) 
b(2.2}»x(3) 
b(3,2)»X(4) 
b(4,2)-X(5) 

a(l,l)*-(a(l,2)+b(l,l)) 
a(2,2) —(a(2,l)+a(2,3)+b(2,2)) 
a(3,3) —(a(3,2)>a(3,4)+b(3,2)) 
a(4,4) —(a(4,3)+a(4,5)-<-b(4,2)) 
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Set initial conditions 


t=0.0 

do i=l,neq 
y(i)=0.0 

do j=l,61 
ys(i,j)=0.0 
enddo 

enddo 

tol=0.0005 

call sset (maxparm, 0.0, param, 1} 
param(4)=2000 


id0=l 

do istep=l,61 

tend=0.07 68 loat{istep) 

CALL DIVPRK (ido, neq, fen, t, tend, tol, param, y) 

do is»l,5 

ys(i,istep)«y(i) 

enddo 

enddo 

Final call to release workspace 


id0a3 

call divprk (idO,neq,fen,t,tend,tol,param,y) 

do i-l,6i 

f(i)-ys(2,i)-t2(i) 

f(i)-ys(3,i)-*t3(i) 

f(i)*ys(4,i)-t4(i) 

enddo 

ssqr-O.O 
do 1^1,61 

ssqr»&sqr4‘f (i) * f (i) 
enddo 

ssqr»ssqr/6l 

xer»sqrt(S8qr) 

write(6,*) xer 

return 

end 


subroutine fcn(neq,t,y,yprinie} 
integer neq 

realms t,y(neq)^yprime(neq) 

real*8 a(5,5),b(5,5),u(5),d,ys(5,6l) 

CQmmon/datal/a,b,u,t2,t3,t4,ys 

if (t.gt.0.2) then 
d-l.O 

else 


d»0.0 
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do i=l,neq 
yprime(i)=0.0 

do j=l,neq 

ypriine(i)=yprinie(i)+a(i, j)*y(j)+b(i, j)*u(j)*d 
enddo 

enddo 


return 

end 





APPEKIDIX D. CONSTRAINED OPTIMIZER - DBCLSF 


This optimizer is the constrained optimizer which allows 
the user to place constraints on the range and values of the 
unknowns. The constraints may be a set of boundary limits or 
restraining the unknown variables to either being negative or 
positive. 


90 





882 Optimization 


BCLSJ/DBCLSJ (Single/Double precision) 

Purpose: Solve a nonlinear least squares problem subject to bounds 

on the variables using a modified Levenberg-Marquardt 
algorithm and a user-supplied Jacobian. 

Usage: CALL BCLSJ (FCN, JAC, M, M. XGUESS, IBTYPE, XLB. XUB, 

XSCALE, FSCALE, IPARAM, RPARAM, X. FVEC, 

FJAC, LDFJAC) 

Arguments 

FCN - User-supplied SUBROUTINE to eveduate the function to be 
minimized. The usage is 
CAa FCH (M. M. X. F). where 

M - Length of F. (Input) 

N - Length of X. (Input) 

X - The point at which the function is evaluated. 
(Input) 

X should not bo changed by FCN. 

F « The computed function at the point X. 

(Output) 

FC?I must be declared EXTERNAL in the calling program. 

JAC - User-supplied SUBROUTINE to evaluate the Jacobian at a 
point X. The usage is 
CALL JAC (M. H, X. FJAC, LDFJAC). where 
M - Length of F. (Input) 

N - Length of X- (Input) 

X - The point at which the function is evaluated. 
(Input) 

X should not be changed by FCN. 

FJAC - The computed M by N Jacobian at the point X. 
(Output) 

LOFJAC - Leading dimension of FJAC. (Input) 

JAC must be declared EXTERNAL in the calling program. 

N - Number of functions. (Input) 

N - Number of variables. (Input) 

XGUESS - Vector of length N containing the initial guess. (Input) 
IBTVPE - Scalar indicating the types of bounds on variables. 

(Input) 

IBTYPE Action 

0 User will supply all the bounds. 

1 All variables are nonnegttlve. 

2 All variables are nonpositive. 

3 User supplies only the bounds on 1st variable, 
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all other variables will have the same bounds. 

XLB - Vector of length N containing the lower bounds on 

variables. (Input, if IBTYPE * 0; output, if IBTYPE » 1 
or 2; input/output, if IBTYPE ■ 3) 

XUB Vector of length M containing the upper bounds on 

variables. (Input, if IBTYPE ■ 0; output, if IBTYPE » 1 
or 2; input/output, if IBTYPE * 3) 

XSCALE - Vector of length N containing the diagonal scaling matrix 
for the variables. (Input) 

In the absence of other information, set all entries to 

1 . 0 . 

FSCALE - Vector of length M containing the diagonal scaling matrix 
for the functions. (Input) 

In the absence of other information, set all entries to 

1 . 0 . 

IPARAH - Parameter vector of length 6. (Input/Qutput) 

See Remarks. 

RPARAM ~ Parameters vector of length 7. (Input/Output) 

See Remarks. 

X - Vector of length N containing the approximate solution. 
(Output) 

FVEC - Vector of length M containing the residuals at she 
approximate solution. (Output) 

FJAC • H by M matrix containing a finite difference approximate 
Jacobian at the approximate solution. (Output) 

LOFJAC.• Leading dimension of FJAC exactly as specified in the 
dimension stattiaent of the calling program. (Input) 

Remarks 

1. Automatic workspace usage is 

BCLSJ 14«N 2*H * I units, or 

08CLSJ 26*N ♦ 4*M - 2 units. 

Uorkspaee say be explicitly provided, if desired, by use of 
B2LSJ/0B2LSJ. The reference is 

CALL B2LSJ (FCN, JAC. M, H. XCUESS, IBTYPE, XLB. XUB. 

XSCALE. FSCALE. IPARAM, RPARAN. X. FVEC. 

FJAC. LSFJAC. WK. IWK) 

The additional arguments are as follows: 

WX - Work vector of length 12«N *- 2«H - 1. WK contains 
the following information on output; 

The second K locations contain the last step taken. 

The third locations contain the last Gauss^Newton step. 
The fourth M locations contain an estimate of the 
gradient at the solution. 


IMSL MATH/LIBRARY 
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IWK - Work vector of lengtii 2«N containing the 

permutations used in the QR factorization of the Jacobian 
at the solution. 

2. Informational errors 
Type Code 

3 1 Both the actual and predicted relative reductions in the 

function are less than or equal to the relative function 
convergence tolerance. 

4 2 The iterates appear to be converging to a noncritical 

point. 

4 3 Maximum number of iterations exceeded. 

4 4 Maximum number of function evaluations exceeded. 

4 5 Five consecutive steps have been taken with the maximum 

step length. 

4 6 Maximum number of Jacobian evaluations exceeded. 

3. The first stopping criterion for BCLSJ occurs when the norm of the 
function is less than Che absolute function colerancs. The second 
stopping criterion occurs when the norm of the scaled gradient is 
less chan the given gradient tolerance. The third stopping 
criterion for BCLSJ occurs when the scaled distance between the 
last two steps is less than the step tolerance. 

4. If nondcfault parameters are desired for IPAIUM or RPAAAK, then 
U4LSF is called and the corresponding parameters are sec to the 
desired value before calling the optimization program. Otherwise, 
if Che default parameters are desired, then set IPARAHd) to zero 
and call the optimization program omitting the call to U4LSF. 

The call to OdLSF would be as follows: 

CALL U4LSF (IPARAN, aPAAAM). 

The following is a list of the parameters and the default values: 
XPAAAM • Integer vector of length 6. 

IPARAMd) Initialization flag. (0) 

IPAftAN(2) • Mumber of good digits in the function. 

(Nacbiae dependent) 

IPAIUNO) ■ .Maximum number of icaraclons. (100) 

IPAIUN(4} ■ Maximaa number of function evaluations. (400) 
IPAIUN(S) ■ Maximum number of Jacobian evaluations. (100) 
IPABiN(6) *• Internal variable scaling flag. (1) 

If IPARAM(6) ■ I Che values for ISCALE are 
set internally. 

aPARAM ■* Real vector of length 7. 
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RPARAM(l) « Scaled gradient tolerance. (eps*«(2/3)) 
RPAilAM(2) • Scaled step tolerance. (eps*»(2/3)) 

RPARAM(3) ■ Relative function tolerance. 

(MAX(1.OE-10.eps»*(2/3))) 

RPARAM<4) • Absolute function tolerance. 

(MAX(1.OE-10,eps**(2/3))) 

RPARAM(S) » False convergence tolerance. (100*eps) 
RPARAM(6) • Maximum allovable step size. 

(1000-MAX(T0H,T0L2)) where. 

TQLl • SQRKsum of (XSCALE(I)-XGUESS(I))—2) 
for I ■ l,...,M 
TQL2 • 2-norm of XSCALE. 

RPARAM(7) « Size of initial trust region radius. 

(Based on the initial scaled Cauchy stop) 
eps is machine epsilon. 

If double precision is desired, then 0U4LSF is called and RPARAM 
is declared double precision. 

Keywords: Levenberg-Marquardt; Trust region 


Algorithm 

3CLSJ uses a modified Levenberg<Marquardt method and an active set strategy to 
solve nonlinear least squares problems subject to simple bounds on the variables. 
The problem is stated as follows: 

mmjFU)^F{r)=5f;A(x)= 

subject to / < X < u. 

where m > n. F : R** — R"*. and /,{x) is the t-th component function of F(x). 
From a given starting point, an active set lA. which contains the indices of the 
variables at their bounds, is built. variable is called a free variable' if it is not in 
the active set. The routine then computes the search direction for the foee variables 
according to the formula 

da 

where p is the L«veaberg*Marquardt parameter. F » F(x). and J is the Jacobian 
with respect to the free variables. The search direction for the variables in lA is 
set to zero. Ttie trust region approach discussed by Dennis and Schnabel (19d3) 
is used to find the uew point. Finaity. the optimaiity conditions are checked. The 
conditions are: 

Ifoi-C.HI < •!. f. <x, <«, 
yix.i <0. r, = (t, 
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gix^) >0. ar, = li, 

where e is a gradient tolerance. This process is repeated until the optimality criterion 
is achieved. 

The active set is changed only when a free variable hits its bounds during an 
iteration, or the optimality condition is met for the free variables but not for adl 
variables in lA, the active set. In the latter case, a variable which violates the 
optimality condition will be dropped out of lA. For more detail on the Levenberg- 
Marquardt method, see Levenberg (1944). or Marquardt (1963). For more detailed 
information on active set strategy, see Gill and .Murray (1976). 


Example 

The nonlinear least squares problem 


min 



“ »al 


subject to -2 < Xi < 0.5 

-1 <-C2 < 2. 

where /i(x) » I 0 (x 2 - xf). and /jfx) » (I - xi) is solved with an initial guess 
(-1.2.1.0). and default values for parameters. 


C 

c 


c 

C 


c 

c 

c 



Oeclaracloa of variables 

ISTEQER LOFJAC. H. H 
PAAAKETER (LI)PJAC«2. K-2. .4«2) 


iKTEGEa mKwa). itp, sour 

(lEAt FJAC(10FJAC..Y). FSCAUfK). PVCC(M). ROSBCX. iUlSJAC. 

k aPAKAN(7). X(M}. XQUESSd). XLa(8}. IS(it). XUB(.'I) 

eXTEUIAL SCLSJ. ROSSa. aOSJAC. tnUCH 

Cocpute Che least squares for the 
RoseabFoch fuact;oa 


DATA XGUESS/-1.2E0. l.OEO/. XS/2-I.OEO/, FSC' 

DATA XL3/-2.0E0. -1.060/, XtIB/O.SEO. 2.0E0/ 

All the bouads are (>rovide<l 


ITP - 0 


IPABANd) • 0 


Default parsMCers are used 


CAU. BOSJ (IbOSBCS. RfiSJAC. M. !«. {GUESS, ITP, XL3. XUB. IS, 
A FSCALE, IPABAM. RPARAN. X. FVEC. FJAC. LOFJAC) 

Print remits 

CALL UNACU (2, fODT) 

vaiTE (ROUT.09999) X. FVEC. IPARANO). IPA11AM(4) 


99999 -OaiUT ' The solutloa is 2F9.4. //. ' The function '. 


aCLSJ>‘DBCLSJ 
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i ‘evaluated at the solutioa is 18X, 2F9.4, //, 

4 ' The auober of iteratioas is lOX, 13, /, ' The 

t ‘niuBber of function evaluations is 13, /) 

END 

C 

SUBROOTINE ROSBCX (M, K, X. F) 

INTEGER N, N 

REAL X(N), F<M) 

C 

F(l) - 1.0E1-(X(2)-X(1)-X(1)) 

F(2) » l.OEO - X(l) 

RETURN 

END 

C 

SUBROUTINE RQSJAC (M. N. X, FJAC, LOFJAC) 

INTEGER M, N. LOFJAC 

REAL XCN). FJAC(LOFJAC.N) 

C 

FJAC(l.l) • -20.0E0*X(l) 

FJAC(2.1) • -l.OEO 
FJAC(l,2) • lO.OEO 
FJAC(2.2) • O.OEO 
RETURN 
END 


Output 


The solutioa is .5000 .2500 

Th« foaetioa svaluatsd at Um solution is 
.0000 .5000 

Th« auober of iteratioas is 13 

the number of function evaluations is 2i 


References 

Dennis. -I. £.. Jr., and Robert 8. Schnabei (1983). Numerical Methods for L'n- 
constrained Optimit»tioa and SoaUnear Equations. Prentice*Hall. Englewood 
CUi&. Mew Jersey. 

Gill. Philip E.. and Walter Murray (1976). Minimaatioa subject to bounds on the 
variables. MPL Report .VAC 72. National Physical Laboratory. England. 

Levenberg, K. (1944). .A method for the solution of irertain problems in (east squares. 
Quarterly of .'{pplied .Mathematics. 2. 164-168. 

Marquardt. 0. (1963). .An algoritiun for least*.<tquares estimation of nonlinear pa¬ 
rameters. SL\M Journal (ut Applied .Mathematics. II. 431-441. 


LMSL m.ath;ubr.arv 


96 


3CLSJ/DBCUJ 





APPENDIX E. ZXSSQ OPTIMIZER 


This optimizer routine was the initial optimizing program 
in the vane modeling process using the IMSL subroutines. It 
provides no way of constraining the values of the unknown 
variables. 
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I.MSL ROUTINE 
PURPOSE 

USAGE 

ARCUHENTS 



NAME - ZXSSQ ' 

< MINIMUM OP the' sum OP SQUARES OF M FUNCTIONS 
IN N VARIABLES USING A FINITE DIFFERENCE 
LEVENBERG-MARQUAROT ALGORITHM 

- CALL ZXSSQtFUNC,M,N,NSIC,EPS,DELTA,MAXFN,lOPT, 
PARM,X.SSQ.F.XJAC.IXJAC.XJTJ,WORK,INFER.lER) 

FUNC - A USER SUPPLIED SUBROUTINE WHICH CALCULATES 
THE RESIDUAL.VECTOR F(1),F(2F(M) FOR 
GIVEN PARAMETER VALUES X(1).X(2X( N). 
THE CALLING SEQUENCE HAS THE FOLLOWING FORM 
CALL FUNC(X.M,N,F) 

WHERE X IS A VECTOR OF LENGTH N AND F IS 
A VECTOR OP LENGTH M. 

FUNC MUST APPEAR IN AN EXTERNAL STATEMENT 
IN THE CALLING PROGRAM. 

FUNC MUST NOT ALTER THE VALUES OF 
X(I).1-1.N, M, OR N. 

H ^ THE NUMBER OF RESIDUALS OR OBSERVATIONS 

(INPUT) 

N - THE NUMBER OF UNKNOWN PARAMETERS (INPUT). 

NSIQ - FIRST CONVERGENCE CRITERION. (INPUT) 

CONVERGENCE CONDITION SATISFIED IF ON TWO 
SUCCESSIVE ITERATIONS, THE PARAMETER 
ESTIMATES AGREE, COMPONENT BY COMPONENT. 

TO NSIG DIGITS. 

EPS • SECOND CONVERCENCC CRITERION. (INPUT) 

CONVERGENCE CONDITION SATISFIED IF. ON TWO 
SUCCESSIVE ITERATIONS THE RESIDUAL SUM 
OF SQUARES ESTIMATES HAVE RELATIVE 
OIFFCRENCe LEv<;S THAN OR EQUAL TO EPS. CPS 
MAY BE SET TO YCRO. 

DELTA - THIRD CONVERGENCE CRITERION. (INPUT) 

CONVERCENCC CCNDITION SATISFIED IF THE 
(EUCLIDEAN) NORM OF THE APPROXIMATE 
GRADIENT IS LESS THAN OR EQUAL TO DELTA. 
DELTA HAY BE SET TO ZERO. 

NOTE. THE ITERATION IS TERMINATED, AND 
CONVERGENCE IS CONSIDERED ACHZ£VCD« IF 
ANY ONE OF THE THREE CONDITIONS IS 
SATISFIED. 

NAXFN • INPUT MAXIMUM NUMBER OP FUNCTION EVALUATIONS 
(I.E., CALLS TO SUBROUTINE FUNC) ALLOWED. 

THE ACTUAL NUMBER OF CALLS TO FUNC MAY 
EXCEED MAXFN SLIGHTLY. 

lOPT INPUT OPTIONS PARAMETER. 

lOPT-0 IMPLIES BROWN'S ALGORITHM WITHOUT 
STRICT DESCENT IS DESIRED. 
tOPT>l IMPLIES STRICT DESCENT AND DEFAULT 
VALUES FOR INPUT VECTOR FARM ARE DESIRED. 
lOPt-2 IMPLIES STRICT DESCENT IS DESIRED WITH 
USER PARAMETER CHOICES IN INPUT VECTOR FARM. 

PABM ' INPUT VECTOR OF LENGTH N USED ONLY FOR 
lOPT EQUAL TWO. PARMd) CONTAINS, WHEN 

ZEROS AND EXTREMA ZXSSQ-1 
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!•), THE INITIAL VALUE OF THE MARQUAROT 
PARAMETER USED TO SCALE THE DIAGONAL OF 
THE APPROXIMATE HESSIAN MATRIX. XJTJ, 

BY THE FACTOR (1.0 ♦ PARM(l)). A SMALL 
VALUE GIVES A NEWTON STEP, WHILE A LARGE 
VALUE GIVES A STEEPEST DESCENT STEP. 

THE DEFAULT VALUE FOR PARM(l) IS 0.01. 

1-2. THE SCALING FACTOR USED TO MODIFY THE 
MARQUAROT PARAMETER, WHICH IS DECREASED 
BY PARH(2) AFTER AN IMMEDIATELY SUCCESSFUL 
DESCENT DIRECTION. AND INCREASED BY THE 
SQUARE OF PARH(2) IF NOT. PARM(2} MUST 
BE GREATER THAN ONE, AND TWO IS DEFAULT. 
1-3. AN UPPER BOUND FOR INCREASING THE 
MARQUAROT PARAMETER. THE SEARCH FOR A 
DESCENT POINT IS ABANDONED IF PARH(3) IS 
EXCEEDED. PARH(3) GREATER THAN 100.0 IS 
RECOMMENDED. DEFAULT IS 120.0. 

I-H. VALUE FOR INDICATING WHEN CENTRAL 

RATHER THAN FORWARD DIFFERENCING IS TO BE 
USED FOR CALCULATING THE .lACOPIAN. THE 
SWITCH IS MADE WHEN THE NORM OF THE 
GRADIENT OF THE SUM OF SQUARES FUNCTION 
BECOMES SMALLER THAN PARM(V), CENTRAL 
DIFFERENCING IS GOOD IN THE VICINITY 
OF THE SOLUTION. SO PARH(V) SHOULD BE 
SMALL. THE DEFAULT VALUE IS 0.10. 

X - VECTOR OP LENGTH N CONTAINING PARAMETER 

VALUES. 

ON INPUT. X SHOULD CONTAIN THE INITIAL 
ESTIMATE OP THE LOCATION OF THE MlNIHUH. 

ON OUTPUT. X CONTAINS THE FINAL ESTIMATE 
OF THE LOCATION OF THE MINIMUM. 

SSQ - OUTPUT SCALAR WHICH 15 SET TO THE RESIDUAL 
SUMS OF SQUARES. F11)•*2♦...•F(M)**2, FOR 
THE FINAL PARAMETER ESTIMATES. 

F - OUTPUT VECTOR OF LENGTH H CONTAINING THE 

RESIDUALS FOR THE PINAL PARAMETER ESTIMATES. 

XJAC - OUTPUT M BY N MATRIX CONTAINING THE 

APPROXIMATE JACOBIAN AT THE OUTPUT VECTOR X. 

XXJAC - INPUT ROW DIMENSION OF MATRIX XJAC EXACTLY 
AS SPECIFIED IN THE DIMENSION STATEMENT 
IN THE CALLING PROGRAM. 

XJTJ - OUTPUT VECTOR OF LENGTH (N*l)*N/2 CONTAINING 

THE N BY N MATRIX (XJAC-TRANSPOSCD) • (XJAC) 
IN SYMMETRIC STORAGE MODE. 

WORK - WORK VECTOR OF LENGTH 5«N ♦ 2*M ♦ {N*l)»N/2. 

ON OUTPUT. WORK(I) CONTAINS FOR 

1-1. THE NORM OP THE GRADIENT DESCRIBED 
UNDER INPUT PARAMETERS DELTA AND PARH(A). 
1-2, THE NUMBER OP FUNCTION EVALUATIONS 
REQUIRED DURING THE W0RK(5) ITERATIONS. 
1-3. THE ESTIMATED NUMBER OF SIGNIFICANT 
DIGITS IN OUTPUT VECTOR X. 

I-)l. THE FINAL VALUE OF THE MARQUAROT 

SCALING PARAMETER DESCRIBED UNDER PARM(l). 

ZEROS AND EXTREMA 
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1*5. THE MUHBER OF ITERATIONS (I.E.. CHANCES 
TO THE X VECTOR) PERPORHEO. 

SEE PROCRAHHINC NOTES FOR DESCRIPTION OF 
THE LATTER ELEHEMTS OF WORK. 

INFER - AM INTEGER THAT IS SET, ON OUTPUT, TO 

INDICATE WHICH CONVERGENCE CRITERION WAS 
SATISFIED. 

INFER - 0 INDICATES THAT CONVERGENCE FAILED. 
lER GIVES FURTHER EXPLANATION. 

INFER • t INDICATES THAT THE FIRST CRITERION 
WAS SATISFIED. 

INFER • 2 INDICATES THAT THE SECOND CRITERION 
WAS SATISFIED. 

INFER - a INDICATES THAT THE THIRD CRITERION 
WAS SATISFIED. 

IF HORE THAN ONE OF T»E CONVERGENCE CRITERIA 
WERE SATISFIED ON THE FINAL ITERATION, 

INFER CONTAINS THE CORRESPONDING SUM. 

(E.C., INFER • 3 IMPLIES FIRST AND SECOND 
CRITERIA SATISFIED SIMULTANEOUSLY). 

ICR > ERROR PARAMETER (OUTPUT) 

TERMINAL ERROR 

ieR»l29 IMPLIES A SINGULARITY WAS DETECTED 
IN THE JACOBIAN AND RECOVERY FAILED. 
ICR«t30 IMPLIES AT LEAST ONE OF M, N, lOPT, 
PARM(I), OR PARM(2) WAS SPCCIFSCD 
INCORRECTLY. 

lER-nZ IMPLIES THAT AFTER A SUCCESSFUL 
RECOVERY FROM A SINGULAR JACOBIAN. THE 
VECTOR X HAS CYCLED BACK TO THE 
FIRST SINGULARITY. 

ICR*133 IMPLIES THAT MAXFN WAS EXCEEDED. 

WARNING ERROR 

IER«36 IMPLIES THAT THE JACOBIAN IS ZERO. 

THE SOLUTION X IS A STATIONARY POINT. 
ICR«39 IMPLIES THAT THE MARQUARDT 
PARAMETER CXCEEDEO PARM(3). THIS 
USUALLY MEANS THAT THE REQUESTED 
ACCURACY WAS NOT ACHIEVED. 

PRECISION > SINGLE AND DOUBLE 
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ZXSSQ Is a finite difference, Levenberg-Marquardt routine for solving 
nonlinear least squares problens. The problem Is stated as follows: 


given H nonlinear functions f^ 

■Inlmlze 
over X ■ 

where x • (x^,X 2 i>>>»X|^) 
be estimated. 


f 2 .of a vector parameter x. 

Is'a vector of N parameters to 


When fitting a nonlinear model to data, the functions f. should 
be defined ae follows: * 

fl (x) ■ 1-1,2,...,M (l.e., the realtjuals) 

where Is the l^th observation of the dependent variable. 

v^ • <vJ,V 2 .''kv^ ** * vector containing the 1-th obser¬ 

vation of the NV independent variables, and 

g la the function defining the nonlinear model. 

ZXSSQ Is based on a modification of the Levenberg-Harquardt algorithm 
which eliminates the need for explicit derivatives. 

Let x^ be an Initial estimate of x. k sequence of approximations 

e» a» 

to the minimum point is generated by 

. • nnnn n«. 

where J. la the numerical Jacobian matrix evaluated at x*^ 

fl • 

0^ la a diagonal matrix equal to the diagonal of 
for lOfTfiO. 

la a positive scaling constant (Narquardt parameter) 

Ifhen forward differences are used, the Jacobian la calculated by 

J 1/2 

where Uj la the J-th unit vector and hj-max(|xj | ,0.l)*eps 

(sec scaling comments in Programming Note 2), where eps la the relative 
precision of floating point arithmetic. For central differences. 

is used. To minimise the number of function evaluations required 
for nonzero lOPt, a rank one update to the Jacobian matrix la used 
where appropriate: 

ZEIOS SNO EXTRENA 
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V, • •’n * TO 

where < - x*'*' - x”. 

See references: 

1 . Brown, K. H., and Dennis. £•. "Derivative free analogues 

of the Levenberg^Marquardt and Gauss algorithms for nonlinear 

least squares approximations", Numerische Mathematlk, i8, 1972, 
289-297. 

2. Brown, K. N., "Computer oriented methods for fitting tabular 
data in the linear and nonlinear least squares sense". Department 
of Computer, Information, and Control Sciences, TR No. 72-13, 
University of Minnesota. 

3. Levenberg, K., "A method for the solution of certain non-linear 

problems in least squares". Quart . Appl . Hath ., 2, I6il-l68. 

Harquardt, D. U., "An algorithm for least-squares estimation 
of nonlinear parameters", SlAM , 1 1 (2)1963. 

Programming Notes 

1. Accuracy in evaluating the functions f^ Is critical when highly 

accurate parameter values, Xy are required (l.e., when NSIC 

la greater than 3 for single precision). In these oases, it 
la advisable to evaluate the functions in precision greater 
than the working precision (e.g., for single precision ZXSSQ, 
double precision can be used). 

2. In reference to significant digit tests for the x vector and 

the sum of squares, leading zeroes to the right of the decimal 
point are counted. For example, both 123.^5 and 0.001 23 have 
five significant digits. Scaling the x vector in the functions 

f^ may be required if several of the parameters x^ are much 

less than 1.0 to obtain the same number of significant digits 
for each of the Xj. 

3. For lOPT 4 0, strict descent Is enforced by increasing the 

Harquardt parameter (see PARM(I)) by the factor PARM(2)^ 

until a smaller sue of squares la obtained. If a smaller S3Q 
Is obtained immediately during Iteration, la decreased by 

the factor PARM(2). also modified by the ratio of 

the norm of the gradient 2||J f|| for two successive iterations 
In the range (PARH(2)'‘\PARH(2)), decreasing as the gradient 

decreases, and increasing gradient increases. 
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For IOPT-0, a formula la used to calculate o , and no check 

la made to see that SSQ has been reduced. The successive Iterates, 

are accepted even if SSQ increases. IOPT-0 seems to be 

moat effective when SSQ has a minimum near zero. In other 
cases and when IOPT-0 fails. It is advisable to use lOPT-i 
or 2. 


5. The Jacoljlan has typical element 3f /Jx., for 1-1,2,...,M and 

* j 2 

J-1,2,...,N. The Hessian has typical element 3 (SSQ ) / 3Xj 3Xj^, 
for j ,k-l ,2 ,. .. N. liote that SSQ-f^f. For small residuals 
f^, the output matrix XJTJ Is a good approximation to the Hessian 
times one half. 


6 . For normal completion (lER-O), a Jacobian calculated by forward 
or central differences (not a rank one update version) Is returned 
In ZJAC. Also, no convergence teats are performed while the 
rank one approximation la being used. 


7. If a terminal error other than IER -130 is returned, all output 

parameters contain the values Indicated for the Iteration during 

which the terminal error occurred. In particular, the five 

values returned In vector WORK should be inspected whether 

a terminal error was returned or not. If 1ER-I3i la returned, 

W0RK(3) will overeatlmate the number of significant digits 

In X by 2 or 3 digits. 

«» 


a. 

9. 


Automatic recovery from singularities in the Jacobian (l.e., 
an entire column is zero) is provided. 


Theoretically, the norm of the gradient. 2||J r)|, should be 
zero at the least squares solution, so that W0HK(1) should 
be inspected before a solution is accepted. Usually, values 


-it 

less than 10 are acceptable, but this number is also dependent 
on the scale of the functions. The graoient times one half 
is located in output vector U0RK(J*K}, where K-KAXKN^l )*N/2,S) 


and J- 1 ,2,...,N. WORK 




3(SSQ}/3x 


J* 


Example 


Suppose the nonlinear model g(x;v^) - X|*v^ - sinlx^ vM is to be 

fitted to 5 data points y^ • v^ ♦ ain(v‘) • (0.1)*(-1)^ where v*-l, 
for i-l.2,...,S. Then f • y,“g{x,v*) and for 

• 4 * 
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APPENDIX F. DATA PROCESSING PROGRAMS 

The programs in this appendix, TRANS.FOR and PRO.M, are 
used to process the experimental data for use in the modeling 
program. TRANS.FOR is a FORTRAN code program which averages 
the values of the experimental data and PRO.M is a MATLAB 
program which normalizes the resultant data, from TRANS.FOR, 
Co ambient and converts the temperatures from degree 
Fahrenheit to degrees Kelvin. 
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Program TRANS 


c 

c This program translates data into free format and smooths 

c the data with a moving window, 

c 

real*8 x(82),y(82),xl(61),yl(61) 
c 

open(8,name= • pre3.dat ', status= • old') 
open(9,name= • tvl.mat ', status® * new') 
open(7,name® • tv2.mat ', status® • new •) 
open(11,name®'tv3.mat',status® * new') 
open(10,name® ' datam.dat •, status® ' new •) 


c Translate the Tvl data 

do i®l,61 
read(8,*)pl 
xl(i)=pl 
1020 enddo 

yl(l)»xl(l) 
yl(2)®xl(2) 
do i=3,59 

yla)=0.2*(xl(J-2)+xl(i-l)^-xl(i)+xl(i+l)+xl{i+2)) 

enddo 

yl(60)=xl(6C; 

yj(6l)®xl(61) 

do i-1,61 

write(9,I000)yl(i) 
c write( 10 ,’k) yl(i) 

enddo 

1000 foriDat(F10.2) 

c Translate the Tv2 data 

do i-1,61 
read(8,*)pi 
xl(i)“pl 
1030 enddo 

yl(l)-xl(l) 
yi(2)«xl(2) 
do i®3,59 

yl (i)-0.2* (xl (i-2)+X1 (i-l) •►xl (i)+X1 (i+i)+xl: i-t-2)) 

enddo 

yl(60)®xl(60) 

yl(6l)«xl(6l) 

do i®l,6l 

write(7,l0\0)yl(i) 
write(l0,*) yl(i) 

enddo 


1010 
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c Translate the Tv3 data 


do i=l,61 
read(8,*)pi 
xl(i)=pl 
1040 enddo 

yl(l)»xl(l) 
yl(2)=xl(2) 
do i=3,59 

y1(i)=0.2*(xl(i-2)+xl(i-1)+xl(i)+xl(i+l)+xl(i+2)) 

enddo 

yl(60)=xl(60) 

yl(61)=xl(61) 

do i»l,61 

write(11,1130)yl(i) 
write(10,*) yl(i) 

enddo 

1130 fonnat(F10.2) 


close(8) 
close(9) 
close(lO) 
close(11) 

end 
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% This program PRO.M normalizes the smoothed data in TV.MAT and TM.MAT 
% to the ambient value which is the first temp in the series of each 
% and then converts the data from degrees F to degrees K. 

% NOTE - THIS PROGRAM IS RUN AFTER THE TRANS.FOR PROGRAM IS RUN 

load tvl 
load tv2 
load tv3 
% 

tvsl=tvl-tvl (1); 
tvsl=tvsl*0. 535; 

% 

tvs2=tv2-tv2(1); 
tvs2=tvs2*0. 555; 

% 

tvs3=tv3-tv3(1); 
tvs3=tvs3*0. 555; 

% 

datam= [ tvsl • tvs2 • tvs3 
datam 

save datam.dat datam /ascii 
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